Spatial constraint based triangular mesh operations in three dimensions

ABSTRACT

A method and system provide the ability to design a (land) surface. A triangular surface mesh representative of an existing surface is obtained. The mesh includes triangles that are connected by vertices and edges. Design constraint sets are determined based design constraints. The design constraints include a maximum slope constraint for a first triangle of the two or more triangles in the triangular surface mesh. The maximum slope constraint is a maximum angle between a normal vector of the first triangle and a reference vector. Heights of the vertices of the first triangle are projected onto the design constraint sets such that the normal vector satisfies all of the design constraints. The projecting includes modifying the heights by a minimum Euclidian distance. A design of the surface represented by the triangular surface mesh is generated based on the projecting.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit under 35 U.S.C. Section 119(e) ofthe following commonly-assigned U.S. provisional patent application(s),which is/are incorporated by reference herein:

Provisional Application Ser. No. 62/394,608, filed on Sep. 14, 2016,with inventor(s) Hung M. Phan and Valentin R. Koch, entitled “ApplyingGeometric Constraints on Triangles in Three Dimensions,”.

BACKGROUND OF THE INVENTION 1. Field of the Invention

The present invention relates generally to road design, and inparticular, to a method, apparatus, system, and article of manufacturefor designing a grading of a triangular mesh representing a physicalsite to accommodate proper drainage.

2. Description of the Related Art

(Note: This application references a number of different publications asindicated throughout the specification by reference numbers enclosed inbrackets, e.g., [x]. A list of these different publications orderedaccording to these reference numbers can be found below in the sectionentitled “References.” Each of these publications is incorporated byreference herein.)

A common representation of three dimensional objects in computerapplications such as graphics and design, are triangular meshes. In manyinstances, individual or groups of triangles need to satisfy spatialconstraints that are imposed either by observation from the real world,or by concrete design specifications for the objects. As these problemstend to be of large scale, choosing a mathematical optimization approachcan be particularly challenging. To better understand the problems ofthe prior art, an explanation of the design problems may be useful.

One issue of the prior art relates to how to solve design problems withconcrete applications in the computational surface generation oftriangular meshes. In many cases, these surfaces need to satisfy certainspatial constraints.

Digital Terrain Modeling in Computational Geometry

Digital terrain modeling may be done through a set of contours, or acollection of ordered elevation points at arbitrary datums in alandscape. These points can then be connected through a triangular mesh,which is usually called a Triangular Irregular Network (TIN) [22]. Aspointed out in [5], traditional triangulation methods may result interrain models that do not represent the true conditions in terms ofdrainage pools and channels, and may result in local minima and sharpslope changes that do not exist in reality. In [6], the authors applyheuristics to generate higher order Delaunay triangulations that attemptto mitigate these problems. A more deterministic approach would be toidentify existing drainage channels and, within certain boundaries,adjust the triangle surfaces in order to smooth out grade change andcreate drainage channels.

Computer-Aided Design of Golf Course Features

Contour changes are an integral part of any golf course design. They mayinvolve smaller terrain changes such as cut and fill of bunker walls forerosion control [13], or higher earthwork quantity features such asartificially embanked mounds, as well as differently sloped greens thatrequire cut and fill in order to fit into the existing ground [13]. Teeflattening with constrained end slopes may require large fillquantities, or can be designed by balancing quantities through benchedfills [15]. Hence, for each feature of a golf course, the architectneeds to consider a variety of surface maximum and minimum slopeconstraints, surface alignment constraints, as well as other surfaceslope requirements that are not only part of the course difficulty, butalso essential for drainage and environmental impact.

Mesh Conversion

A common problem in topology is the conversion of one surfacerepresentation into another. One such example is converting an existingtriangular, irregular mesh into the closest regular rectangular mesh. Ineach context, the term closest could mean different things such as thevertical distance of the corner points of the rectangle to thetriangular grid, or the minimum distance of the average of the triangleelevations to the average of the corner points, or the minimum of thevolume between the triangles and the rectangular surface thatencompasses them. Methods to convert two dimensional, triangular meshesinto quadrilateral meshes while keeping some or all of the originalvertices are shown in [17, 18, 24]. An extension of these methods tothree dimensions would require that the triangles that are grouped intoquadrilaterals to satisfy some surface alignment conditions.Furthermore, one would like to minimize the distance between theoriginal and final elevations of the vertices, using either the 1-normor the 2-norm. FIG. 1 shows an example of an existing triangular,irregular mesh at the bottom, and a quadrilateral grid on the top. Thedesigner would like to convert the three dimensional triangular mesh 102into a three dimensional quadrilateral mesh that aligns with the gridlines 104 shown on the top. FIG. 2 shows how the existing triangularmesh 102 is mapped in the quadrilateral grid 104. If a quadrilateralcuts an existing triangle, it may need to be retriangulated to form newtriangles so that each quadrilateral is triangulated. The mentionedsurface alignment conditions would then apply to triangles that areinside an individual quadrilateral. The design objective function isthen to minimize the elevation difference between the triangle verticesin the original triangular mesh 102, and the triangle vertices in thenew mesh that groups the triangles into quadrilaterals.

Computer-Aided Design for Architecture and Civil Engineering Structures

In Computer-Aided Design (CAD), triangular meshes are widely used invarious engineering disciplines. Existing and finished ground surfacesin architecture and civil engineering are represented by triangulatedirregular networks. Slopes are relevant in the context of drainage, inboth, civil engineering (transportation structures), and architecture(roof designs), as well as in the context of surface alignments such assolar farms, embankment dams, and airport infrastructure layouts.

In architecture, the design of a complex roof requires slope constraintsfor drainage, but can also include orientation and maximum slopeconstraints. Consider the example in FIG. 3 of a roof seen from above.The roof polygons 302 with the light dashed lines are triangulated, inorder to create a triangular mesh. The engineer may want to placemultiple solar panels 304 on a roof. The solar panels 304 require aminimum inclination towards the position of the sun at noon, shown bythe arrows 306. Also, the panel 304 should not exceed a certain angle,to preserve a maximum average efficiency. This requires a surfaceorientation constraint for some of the triangles.

A concrete problem that arises in civil engineering design is thegrading of a parking lot. Within a given area, the engineer has todefine grading slopes such that the parking lot fits within existingstructures, the drainage requirements on the lot are met, and safety andcomfort is taken into account. Besides these requirements, the engineerwould like to change the existing surface as little as possible, inorder to save on earthwork costs.

FIG. 4 depicts an example of a parking lot that has four horizontalparking rows. In FIG. 4, the designer wants water to flow in thedirection of the arrows 402 towards and along the thick lines 404 intothe inlets 406 that are represented as circles. Lines 408 are curbs thatneed to be aligned.

Drainage schemes can vary for each situation. The triangular mesh inFIG. 5 represents the existing ground of a planned parking lot. Theengineer would like the water to drain away from the building, and alongthe drainage lines 502 into the four corners. Lines 504 indicate wherethe triangle edges need to be aligned. In FIG. 6, the engineer wouldlike to study a different scheme, where the drainage 602 happens inparallel, on either side of the building. Lastly, FIG. 7 represents themesh of an existing ground for a roundabout, where a minimum slope isrequired from the inner circle to the outer one, and water needs todrain from the road at the top along the outer circle to the roads oneither side at the bottom (indicated by drainage lines 702).

SUMMARY OF THE INVENTION

Embodiments of the invention overcome the problems of the prior artbased on geometric design constraints for triangles. In particular, adesign constraint is defined that includes a maximum slope constraintfor at least one of the triangles in a triangular surface mesh. Themaximum slope constraint is a maximum angle between a normal vector ofthe triangle and a reference (e.g., Z) vector. Heights of the verticesof the triangle are projected onto design constraint sets such that thenormal vector satisfies all of the design constraints (e.g., the maximumslope constraint). During the projection, the height(s) (e.g., theZ-coordinate [and not the x or y coordinate(s)]) of the triangle ismoved by a minimum Euclidian distance. The projection is performediteratively such that the triangular mesh converges to an optimalsolution for a CAD design of the land surface. Accordingly, a design ofthe surface represented by the triangular surface mesh is generatedbased on the projecting.

BRIEF DESCRIPTION OF THE DRAWINGS

Referring now to the drawings in which like reference numbers representcorresponding parts throughout:

FIG. 1 shows an example of an existing triangular, irregular mesh at thebottom, and a quadrilateral grid on the top;

FIG. 2 shows how an existing triangular mesh is mapped in aquadrilateral grid;

FIG. 3 illustrates an exemplary roof structure view from above;

FIG. 4 depicts an example of a parking lot that has four horizontalparking rows;

FIG. 5 illustrates a triangle mesh representative of an existing groundof a planned parking lot;

FIG. 6 illustrates a triangular mesh for a planned parking lot wheredrainage happens in parallel, on either side of the building;

FIG. 7 represents the mesh of an existing ground for a roundabout, wherewater needs to drain from the road at the top along the outer circle tothe roads on either side at the bottom;

FIG. 8 shows a triangular mesh in accordance with one or moreembodiments of the invention;

FIG. 9 illustrates a simple triangular mesh with vertices, edges, andtwo triangles in accordance with one or more embodiments of theinvention;

FIG. 10 illustrates two lines A and B that are shown as examples of twoconstraint sets in accordance with one or more embodiments of theinvention;

FIG. 11 shows the iterates of the DR method on two constraint sets A andB in accordance with one or more embodiments of the invention;

FIG. 12 shows a triangular mesh from a birds-eye view in accordance withone or more embodiments of the invention;

FIG. 13 shows the birds-eye view of a feature line on a triangular meshin accordance with one or more embodiments of the invention;

FIG. 14 illustrates a least squares problem in accordance with one ormore embodiments of the invention;

FIG. 15 illustrates a maximum angle for a normal vector of a groundplane with respect to gravity in accordance with one or more embodimentsof the invention;

FIG. 16 illustrates an optimal solution based on tangent points betweena circle and level sets in accordance with one or more embodiments ofthe invention;

FIG. 17 illustrates a desired orientation ({right arrow over (q)}) of anormal vector ({right arrow over (n)}) for a triangle plane inaccordance with one or more embodiments of the invention;

FIG. 18 illustrates a visualization of a surface minimum slopeconstraint in accordance with one or more embodiments of the invention;

FIG. 19 illustrates a feasible set region for a normal vector withrespect to a Surface Oriented Minimum Slope Constraint in accordancewith one or more embodiments of the invention;

FIG. 20 shows that there are exactly two tangent points (u, v) withpositive Lagrange multipliers λ, one of which has u<0 and the other hasu>0 in accordance with one or more embodiments of the invention;

FIG. 21 shows that there is exactly one tangent point (u, v) with apositive Lagrange multiplier λ and u>0, and at least one tangent pointwith a nonpositive Lagrange multiplier in accordance with one or moreembodiments of the invention;

FIG. 22 illustrates a circular segment defining a feasible set for amaximum slope and an oriented minimum slope in accordance with one ormore embodiments of the invention;

FIG. 23 illustrates a circle representing a feasible region for amaximum slope constraint and for two oriented minimum slope constraintsin accordance with one or more embodiments of the invention;

FIG. 24 illustrates multiple surface oriented minimum slope constraintsvia the intersection of two corresponding circular segments inaccordance with one or more embodiments of the invention;

FIG. 25 illustrates two adjacent triangles with surface normal vectorsthat are used to minimize the curvature in accordance with one or moreembodiments of the invention;

FIG. 26 illustrates a hexagon that is used to replace a unit disk for amaximum slope constraint in accordance with one or more embodiments ofthe invention;

FIG. 27 shows how the original half cone is now replaced with theabsolute value function in accordance with one or more embodiments ofthe invention;

FIG. 28 shows the starting situation for the problem shown in FIG. 7 inaccordance with one or more embodiments of the invention;

FIG. 29 shows the mesh of FIG. 28 after 10 iterations in accordance withone or more embodiments of the invention;

FIG. 30 shows the mesh of FIG. 28 after 29 iterations in accordance withone or more embodiments of the invention;

FIG. 31 shows the final mesh of FIG. 28 after convergence in accordancewith one or more embodiments of the invention;

FIG. 32 illustrates the mesh of FIG. 28 based on the problem shown inFIG. 8 after 15 iterations in accordance with one or more embodiments ofthe invention;

FIG. 33 illustrates the mesh of FIG. 28 based on the problem shown inFIG. 8 after 100 iterations;

FIG. 34 illustrates the final mesh of FIG. 28 based on the problem shownin FIG. 8 in accordance with one or more embodiments of the invention;

FIG. 35 is the starting situation mesh for the problem shown in FIG. 9in accordance with one or more embodiments of the invention;

FIG. 36 illustrates the mesh for the problem shown in FIG. 9 after 10iterations in accordance with one or more embodiments of the invention;

FIG. 37 illustrates the mesh for the problem shown in FIG. 9 after 20iterations in accordance with one or more embodiments of the invention;

FIG. 38 is the final solution mesh for the problem shown in FIG. 9 after20 iterations in accordance with one or more embodiments of theinvention;

FIG. 39 illustrates the logical flow for designing a surface inaccordance with one or more embodiments of the invention;

FIG. 40 is an exemplary hardware and software environment used toimplement one or more embodiments of the invention; and

FIG. 41 schematically illustrates a typical distributed/cloud-basedcomputer system using a network to connect client computers to servercomputers in accordance with one or more embodiments of the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

In the following description, reference is made to the accompanyingdrawings which form a part hereof, and which is shown, by way ofillustration, several embodiments of the present invention. It isunderstood that other embodiments may be utilized and structural changesmay be made without departing from the scope of the present invention.

The figures herein are used to illustrate various geometric constraintson the triangles, design applications, and the way how iterativeprojections are performed.

Mathematical notation is utilized herein for a precise description ofthe geometric constraints and the corresponding projections. Themathematical notation is fairly standard and follows largely [2]. Rdenotes the set of real numbers. By “x:=y”, or equivalently “y=: x”, wemean that “x is defined by y”. The assignment operators are denoted by“δ” and “→”. The angle between two vectors {right arrow over (n)} and{right arrow over (q)} is denoted by ∠({right arrow over (n)},{rightarrow over (q)}). The cross product of {right arrow over (n)}₁ and{right arrow over (n)}₂ in R³ is denoted by {right arrow over(n)}₁×{right arrow over (n)}₂.

General Overview

Embodiments of the invention provide a method and system for utilizingprojections to various geometric design constraints for triangles inthree dimensions, and proximity operators for design objectivefunctions. The projections and proximity operators are executed ontriangles that are part of an existing triangular mesh. Further, theseprojections and proximity operators are executed in an iterative way,such that the triangular mesh converges to an optimal solution for acomputer-aided design in three dimensions.

Problem Formulation

A common representation of three dimensional objects in computerapplications such as graphics and design, are triangular meshes. In manyinstances, individual or groups of triangles need to satisfy spatialconstraints that are imposed either by observation from the real world,or by concrete design specifications for the objects. As these problemstend to be of large scale, choosing a mathematical optimization approachcan be particularly challenging. For embodiments of the invention,various geometric constraints are defined as convex sets in Euclideanspaces, and the corresponding projections are found. Furthermore,proximity operators are introduced for certain objective functions. Theprojections and the proximity operators are used in an iterative methodto find optimal design solutions.

Triangular Mesh

In computer-aided design, a triangular mesh is a surface that isrepresented by a set of connected triangles in three dimensional spaceR³. This form of mesh is also called a Triangular Irregular Network(TIN). FIG. 8 shows a triangular mesh in R³ with underlying contours onthe R² plane in accordance with one or more embodiments of theinvention. What follows is a precise mathematical definition of such atriangular mesh.

Assume that G=(V₀,E₀) is a given triangular mesh on R² where V₀ is theset of vertices and E₀ is the set of (undirected) edges, i.e.,V ₀ :={p _(i)=(p _(i1) ,p _(i2))∈R ² |i∈{1,2, . . . ,n}},  (1)E ₀ ⊂{p _(i) p _(j) ≡p _(j) p _(i) |p _(i) ,p _(j) ∈V ₀}.  (2)From G, one can derive the set of triangles of the mesh as followsT ₀:={Δ=(p _(i) p _(j) p _(k))p|{p _(i) p _(j) ,p _(j) p _(k) ,p _(k) p_(i) }⊂E ₀}.  (3)

One first aims tofind a vector z=(z ₁ , . . . ,z _(n))∈X  (4)such that the points{P _(i)=(p _(i1) ,p _(i2) ,z _(i))}_(i∈{1, . . . ,n})satisfy a given setof constraints.  (5)

Clearly, the points (P_(i))_(i∈{1, . . . ,n}) also form a correspondingtriangular mesh S in three dimensions, where the mesh does not containtriangles that have a vertical surface (since we first defined all thevertices as unique points in R²). FIG. 9 illustrates a simple triangularmesh with verticesV ₀ :={P ₁ ,P ₂ ,P ₃ ,P ₄},edgesE ₀ :={P ₁ P ₄ ,P ₁ P ₂ ,P ₂ P ₄ ,P ₃ P ₄ ,P ₂ P ₄},and the two trianglesT ₀:={(P ₁ P ₂ P ₄),(P ₂ P ₃ P ₄)}.in accordance with one or more embodiments of the invention.

Therefore, if there is no confusion, E₀, V₀, T₀ will be reused to denotethe vertices, the edges, and the triangles of the three dimensionalmesh.

Interval Constraint

A designer may require that some of the vertices in the triangular meshcannot exceed a certain elevation, or cannot be lower than a certainelevation. Hence, the elevation of these vertices needs to lie within agiven interval. This is referred to as the interval constraint and isdefined as follows. For a given subset I of {1, 2, . . . , n}, for alli∈I, the entries z_(i) must lie in a given interval of R.

Edge Slope Constraint

A designer may require that the slope of a triangle edge cannot besteeper than a given value, or must be steeper than a given value.Hence, the edge slope needs to lie within a given interval. This isreferred to as the edge slope constraint and is defined as follows. Fora given subset E of the edges E₀, and for every edge p_(i)p_(j)∈E, theslope

$\begin{matrix}{s_{ij}:=\frac{z_{i} - z_{j}}{\sqrt{\left( {p_{i\; 1} - p_{j\; 1}} \right)^{2} + \left( {p_{i\; 2} - p_{j\; 2}} \right)^{2}}}} & (6)\end{matrix}$must lie in a given subset of R.

Edge Alignment Constraint

A designer may require that the edges of some of the triangles in themesh align themselves. This is referred to as the edge alignmentconstraint and is defined as follows. For a given pair of edgesp_(i)p_(j), p_(m)p_(n)∈E₀, the edges must have the same slopes_(ij)=s_(mn).

Surface Alignment Constraint

A designer may require that the surfaces of some of the triangles in themesh align themselves. Aligning triangle surfaces can result in usefulpatches inside the mesh that may be non-triangular. This is referred toas the surface alignment constraint and is defined as follows. For agiven pair of triangles Δ_(i), Δ_(j)∈T₀ the normal vectors {right arrowover (n)}_(Δ) _(i) and {right arrow over (n)}_(Δ) _(j) must be parallel.

Surface Slope Constraint

A designer may require that one or multiple triangles cannot exceed acertain slope with respect to gravity. This is referred to as a maximumsurface slope constraint. A designer may also require that sometriangles need to have a certain minimum slope in the direction of agiven vector. This is referred to as an oriented minimum surface slopeconstraint. Both surface slope constraints can be defined as a moregeneral surface slope constraints that follow. For a given subset T ofthe triangles T₀ and for each triangle Δ∈T, there is a given set ofvectors Q_(Δ) ⊂R³ such that for each {right arrow over (q)}∈Q_(Δ), theangle between the normal vector {right arrow over (n)}_(Δ) and {rightarrow over (q)} must lie in a given subset θ_({right arrow over (q)}) of[0, π], which is also referred to as a surface orientation constraint.

Curvature Minimization

In some design problems, the designer may wish to construct a surfacethat is as “smooth” as possible. This problem is referred to asminimizing the curvature between adjacent triangles in the mesh. One wayto minimize the curvature is to minimize the grade change between eachpair of adjacent triangles in the mesh. One may define this as follows.For a given pair of triangles Δ_(i), Δ_(j)∈T₀ the normal vectors {rightarrow over (n)}_(Δ) _(i) and {right arrow over (n)}_(Δ) _(j) we want tominimize the angle between the two normal vectors.

Methods Overview

Embodiments of the invention are based on projections, and the use ofprojections in iterative ways to find solutions for feasibility oroptimization problems. A projection of a point onto a convex set is theshortest distance of that point to the set. For example, if a set is aplane in three dimensions, then the projection of a point onto that setis its shadow on the plane. If the point to be projected is already onthe plane, than the projection is the point itself.

Given multiple sets (e.g. multiple planes in three dimensions) and astarting point, iterative projection methods can be used to find a pointin the intersection of all these sets that is closest from the startingpoint. These methods project iteratively onto the different sets. If thegiven sets are convex, then the iterative projections will converge to asolution.

In embodiments of the invention, the sets for each constraint areidentified. Closed-form projections are then developed for each of theconstraint sets. A proximity operator can be found for the curvatureminimization. Iterative projection methods, such as the Douglas-Rachfordsplitting method, are then used to project iteratively to the constraintsets. This alters the elevations of the triangles in the triangular meshuntil the mesh converges to a configuration that is desired by theengineer.

Projections onto Constraint Sets

Suppose there are J constraints imposed on a model. For j ∈{1, 2, . . ., J}, we denote by C_(j) the set of all points that satisfies the j-thconstraint. Thus, (4) can be written in the mathematical form

$\begin{matrix}{{{{{find}\mspace{14mu} a\mspace{14mu}{point}\mspace{14mu} x} \in C}:={\bigcap\limits_{j \in {\{{1,2,\;\ldots\;,J}\}}}C_{j}}},} & (7)\end{matrix}$which is referred to herein as a feasibility problem. Moreover, of theinfinitude of possible solutions for (7), one may be particularlyinterested in those that are optimal in some sense. For instance, itcould be desirable to find a solution that minimizes the slope changebetween adjacent triangles, a solution that minimizes the volume betweenthe initial triangles and the triangles in the final solution, orvariants and combinations thereof.

If more than one objective function is of interest, it is common toadditively combine these functions, perhaps by scaling the functionsbased on their different levels of importance. In summary, one goal isto solve the problem

$\begin{matrix}{{{{\min\mspace{14mu}{F(z)}\mspace{14mu}{subject}\mspace{14mu}{to}\mspace{14mu} z} \in C}:={\bigcap\limits_{j \in {\{{1,2,\;\ldots\;,J}\}}}C_{j}}},} & (8)\end{matrix}$where F may be a sum of (scaled) objective functions. The function F istypically nonsmooth which prevents the use of standard optimizationmethods.

A constraint set C⊂X is the collection of all feasible data points,i.e., points that satisfy some requirements. Suppose the given datapoint z∈X is not feasible (i.e., z∉C), one aims to modify z so that thenewly obtained point x is feasible (i.e., x∈C ); and one would like todo it with minimal adjustment on z. This task can be achieved by usingthe projection onto C. Recall that the projection of z onto C, denotedby P_(C)z, is the solution of the optimization problem

$\begin{matrix}{{P_{C}z}:={{\underset{\;{x \in C}}{\arg\;\min}{{{z - x}}}} = {\left\{ \;{{x \in C},{{{{x - z}}} = {\min\limits_{\;{y \in C}}{{{z - y}}}}}} \right\}.}}} & (9)\end{matrix}$

It is well known that if C is nonempty, closed, and convex¹, then P_(C)zis singleton, see, e.g., [26, Theorem 2.10]. ¹C is convex if for all x,y∈C and λ∈[0,1], we have (1−λ)x+λy∈C

It turns out that all spatial constraints encountered in variousexemplary settings may be convex and closed. Hence, their projectionsalways exist and are unique. Moreover, one can also successfully obtainexplicit formulas for these projections.

Iterative Projection Methods for Feasibility Problems

Consider the feasibility problem

$\begin{matrix}{{{{{find}\mspace{14mu} a\mspace{14mu}{point}\mspace{14mu} x} \in C}:={\bigcap\limits_{j \in {\{{1,2,\;\ldots\;,J}\}}}C_{j}}},} & (10)\end{matrix}$

The intersection C is the feasible set containing the desired solutions.Clearly, it might be tempting to solve this problem by finding theprojection onto C directly. However, this is often not possible due tothe complexity of C. A workaround is to utilize the projection P_(C)_(j) onto each constraint, if the explicit formula is available. Thenwith an initial point, one can iteratively execute the projections P_(C)_(j) 's to derive a solution for (10). Let z₀ ∈X be the initial datapoint. The following are just a few of many (iterative) projectionmethods (see [2, 10, 12] and the references therein).

Cyclic Projections

In this method, a point is projected in sequence, from one constraintset to the next. Once the projected point reaches the first constraintset again, it is referred to as one cycle or one iteration. FIG. 10illustrates two lines A and B that are shown as examples of twoconstraint sets. The method starts with point b⁻¹ that is projected tothe set A, resulting in the new point a₀ that then gets projected to theset B. The obtained point b₀ completes one cycle. The process isrepeated until the iterations converge to the intersection 1002 of thetwo sets A and B. The cyclic projections are defiend as follows. Givenz_(k), one updates z_(k+1):=Tz_(k) whereT:=P _(C) _(J) P _(C) _(J-1) . . . P _(C) ₂ P _(C) ₁ .  (11)

Parallel Projections

In this method, a point is projected in parallel to all the constraintsets at once. The resulting points are then averaged to obtain a newpoint for the next iteration. One defines the parallel projections asfollows. Given z_(k), one updates z_(k+1):=Tz_(k) where

$\begin{matrix}{T:={\frac{1}{J}{\left( {P_{C_{1}} + P_{C_{2}} + \ldots + P_{C_{J}}} \right).}}} & (12)\end{matrix}$

String-Average Projections

The string-average projections is a combination of cyclic and parallelprojections. One defines it as follows. Given z_(k), one updatesz_(k+1):=Tz_(k) where

$\begin{matrix}{T:={\frac{1}{J}{\left( {P_{C_{1}} + {P_{C_{2}}P_{C_{1}}} + \ldots + {P_{C_{J}}P_{C_{J - 1}}\mspace{14mu}\ldots\mspace{14mu} P_{C_{2}}P_{C_{1}}}} \right).}}} & (13)\end{matrix}$

In fact, these are probably the simplest methods for solving feasibilityproblems. In addition, when projection methods succeed (see [11] forsome interesting examples), they have various attractive features: theyare easy to understand, simple for implementation and maintenance, andsometime can be very fast. [1, 3, 10, 12] provide additional discussionson projection methods.

Iterative Projection Methods for Optimization Problems

Iterative optimization methods are often used for solving (8), which mayrequire the computations of proximity operators for the functionsinvolved. Suppose ƒ:X→(−∞,+∞] is a proper convex lower semicontinuousfunction (see [25] and [2] for a discussion of convex analysis) and fixx∈X. Then it is well known (see [2, Section 12.4]) that the function

$\begin{matrix}\left. X\rightarrow{\left( {{- \infty},{+ \infty}} \right\rbrack\text{:}\mspace{14mu} y}\mapsto{{f(y)} + {\frac{1}{2}{{{x - y}}}^{2}}} \right. & (14)\end{matrix}$has a unique minimizer, which will be denoted by P_(ƒ) (x). The inducedoperator P_(ƒ): X→X is called the proximal mapping or proximity operatorof ƒ (see [23]). Note that if ƒ is the indicator function of a set C(the indicator function t_(C) is defined by t_(C)(x)=0 if x∈C andt_(C)(x)=+∞ otherwise), then P_(ƒ)=P_(C). Thus, proximity operators aregeneralizations of projections.

The Douglas-Rachford Method

As proximity operators may be made available for several types ofobjective functions, any iterative optimization methods that utilizeproximity operators (e.g., [7, 8, 21]) can be used to solve thecorresponding problems. Let us describe one such method, theDouglas-Rachford (DR) splitting. DR emerged from the field ofdifferential equations [14], and was later made widely applicable inother areas thanks to the seminal work [21]. FIG. 11 shows the iteratesof the DR method on two constraint sets A and B in accordance with oneor more embodiments of the invention.

Using indicator functions, (8) is converted intomin F(x)+t _(C) ₁ +t _(C) ₂ + . . . +t _(C) _(J) subject to x∈X.  (15)

So, it suffices to present DR for the following general problem

$\begin{matrix}{{{\min\mspace{14mu}{\sum\limits_{i = 1}^{m}\;{{f_{i}(x)}\mspace{14mu}{subject}\mspace{14mu}{to}\mspace{14mu} x}}} \in X},} & (16)\end{matrix}$where each ƒ_(i) is a proper convex lower semicontinuous function on X.The DR operates in the product space X:=X^(m) with inner product

x,y

:=Σ_(i=1) ^(m)

x_(i), y_(i)

for x=(x_(i))_(i=1) ^(m) and y=(y_(i))_(i=1) ^(m). Set the startingpoint x₀=(z, . . . , z)∈X, where z∈X. Given x_(k)=(x_(k,l), . . . ,x_(k,m))∈X, one computes

$\begin{matrix}{{{\overset{\_}{x}}_{k}:={\frac{1}{m}{\sum\limits_{i \in K}\; x_{k,i}}}},} & (17)\end{matrix}$∀i∈{1, . . . ,m}:x _(k+1,i) :=x _(k,i) −x _(k) +P _(γƒ) _(i) (2 x _(k,i)−x _(k,i)),  (18)then update x _(k+1):=(x _(k+1,l) , . . . ,x _(k+1,m)).  (19)Then the monitored sequence (x _(k))_(k∈N) converges to a solution of(16).

It is worth mentioning that when all ƒ_(i)'s are indicator functions ofthe constraints (thus, all proximity operators become projections), then(16) becomes a pure feasibility problem. Therefore, all optimizationmethods for (16) can also be used for feasibility problems.

Parallelization

By using the indicator function, projections become a special case ofproximity operators. Thus, the following remarks simply refer toproximity operators.

Remark 1 (Entries that Possibly Change after Executing a ProximityOperator)

Suppose X=R^(n) and the objective function ƒ: X→(−∞, +∞] only depends oncertain entries, i.e., ƒ(x)=ƒ((x_(i))_(i∈I)), where I⊂{1, 2, . . . , n}is a given index set.

Let z=(z₁, . . . , z_(n))∈X and execute the proximity operator P_(ƒ)Z.It is obvious that after the execution, the coordinates that possiblychange are (z_(i))_(i∈I); while all other coordinates (z_(i))_(i∉I) areleft unchanged, see, e.g., Method component 2 (projection onto an edgeminimum slope constraint). Thus, for simplicity of the presentation, foreach projection or proximity operator, one may often focus on the“possibly affected” coordinates (z_(i))_(i∈I) in the subspace R^(|I|)where |I| is the number of indices in I.

Remark 2 (Parallelization)

If two functions ƒ₁, ƒ₂: X→(−∞,+∞] are independent in the sense thatthey depend on separate coordinates, then one can execute theirrespective proximity operators in parallel. Therefore, one can packmultiple independent proximity operators to accelerate the convergenceof the method.

Projections onto Linear Constraints

The following describes the projections found for the linear constraintsthat were defined in the problem formulation above. These closed-formprojections ensure that the method will converge to a solution in areasonable amount of time, using an iterative projection method.

By linear constraints, one refers to any constraint on the triangularmesh that can be represented by a system of linear equalities andinequalities. Indeed, this class includes several important constraintsin design problems. In the sections that follow, those constraints willbe analyzed.

Interval Constraint

Similar to [3, Section 2.2], one may assume that I is a subset of {1, 2,. . . , n} and (l_(i))_(i∈I), (u_(i))_(i∈I)∈R^(I) are given. SetY:={z=(z ₁ ,z ₂ , . . . ,z _(n))∈X|∀i∈I:l _(i) ≤z _(i) ≤u _(i)}.  (20)

The closed set Y is an affine subspace, thus, convex. The followingexplicit formula is for the projection onto Y, whose proof isstraightforward, thus, omitted.

Method Component 1 (Projection onto Interval Constraints)

$\begin{matrix}{\left. {P_{Y}\text{:}\mspace{14mu} X}\rightarrow{X\text{:}\mspace{14mu}\left( {z_{1},z_{2},\ldots\mspace{11mu},z_{n}} \right)}\mapsto\left( {x_{1},x_{2},\ldots\mspace{11mu},x_{n}} \right) \right.,} & (21) \\{{{where}\mspace{14mu} x_{i}} = \left\{ {\begin{matrix}{\max\left\{ {l_{i},{\min\left\{ {u_{i},z_{i}} \right\}}} \right\}} & {{i \in I},} \\{z_{i},} & {i \in {\left\{ {1,2,\ldots\mspace{11mu},n} \right\}\backslash I}}\end{matrix}.} \right.} & (22)\end{matrix}$

Edge Minimum Slope Constraint

Let P_(i)=(p_(i1), p_(i2), h_(i)) and P_(j)=(p_(j1), p_(j2), h_(j)). Thedesigner may require that the (directional) slope from P_(j) to P_(i)must be no less than a threshold level s_(ij). More specifically, sincethe value p_(i1), p_(i2), p_(j1), p_(j2) are fixed, one can write thisconstraint ash _(i) −h _(j)≥α_(ij) :=s _(ij)√{square root over ((p _(i1) −p_(j1))²+(p _(i2) −q _(j2))²)},  (23)which is called the edge minimum slope constraint. This constraint is alinear inequality, thus, convex. Projection on to this type ofconstraint is derived analogously to [3, Section 2.3], thus, its proofis also omitted.

Method Component 2 (Projection onto an Edge Minimum Slope Constraint)

LetE _(i,j) ={x=(x ₁ ,x ₂ , . . . ,x _(n))∈R ^(n) |x _(i) −x_(j)≥α_(ij)}  (24)be the minimum slope constraint on the edge P_(i)P_(j). Let x:=(x₁, . .. , x_(n)):=P_(E) _(i,j) z be the projection of a point z=(z₁, z₂, . . ., z_(n)) on to E_(i,j). Then

$\begin{matrix}{{x_{i} = {\frac{1}{2}\left( {z_{i} + z_{j} + {\max\left\{ {{z_{i} - z_{j}},\alpha_{ij}} \right\}}} \right)}},} & (25) \\{{x_{j} = {\frac{1}{2}\left( {z_{i} + z_{j} - {\max\left\{ {{z_{i} - z_{j}},\alpha_{ij}} \right\}}} \right)}},{and}} & (26) \\{x_{k} = {{z_{k}\mspace{14mu}{for}\mspace{14mu}{all}\mspace{14mu} k} \notin {\left\{ {i,j} \right\}.}}} & (27)\end{matrix}$

Low Point Constraint

A point P_(j)=(p_(j1), p_(j2), z_(j)) on the mesh is called a low pointif each point P_(i)=(p_(i1), p_(i2), z_(i)) connected to P_(j) satisfiesthe edge minimum slope constraintz _(i) −z _(j)≥α_(i) :=s _(i)√{square root over ((p _(i1) −p _(j1))²+(p_(i2) −q _(j2))²)},  (28)where all s_(i)∈R are given.

FIG. 12 shows a triangular mesh from a birds-eye view, where the lowpoint P₀ is connected to points P₁, P₂, P₃, P₄, P₅, and P₆. So in thispicture, one requires that P₀ is lower than all its connected points andthat the slopes on each connecting edge i is at least s_(i).

One can treat the low point constraint as a single constraint eventhough it is the intersection of finitely many edge slope constraints.The following result shows how to project onto this constraint.

Method Component 3 (Projection onto a Low Point Constraint)

Let α₂, . . . , α_(m)∈R. Define the setE:={(x ₁ , . . . ,x _(m))∈R ^(m) |∀i∈{2, . . . ,m}:x _(i) −x₁≥α_(i)}.  (29)

Let z=(z₁, . . . , z_(m))∈R^(m). Let δ₁=z₁ and let δ₂≤δ₃≤ . . . ≤δ_(m)be the values {z_(i)−α_(i)}_(i∈{2, . . . , m}) in nondecreasing order.Let k be the largest number in {1, . . . , m} such that δ_(k)≤(δ₁+ . . .+δ_(k))/k. Then the projection x:=(x₁, x₂, . . . , x_(m)):=P_(E)Z isgiven byx ₁=(δ₁+ . . . +δ_(k))/k and x _(i)=max{x ₁ ,z _(i)−α_(i)}+α_(i) ,∀i∈{2,. . . ,m}.   (30)

To show this formula, first, x is the solution ofmin(x ₁ −z ₁)²+ . . . +(x _(m) −z _(m))²  (31)s.t x ₁ −x _(i)+α_(i)≤0,i∈{2, . . . ,m}.  (32)

Let y:=(y₁, . . . , y_(m)) where y₁=x₁ and y_(i):=x_(i)−α_(i) for i∈{2,. . . , m}. Without loss of generality, one may assume δ_(i)=z_(i)−α_(i)for i∈{2, . . . , m}. Then (4) becomesmin ϕ(y):=(y ₁−δ₁)²+ . . . +(y _(m)−δ_(m))²  (33)s.t g _(i)(y):=y ₁ −y _(i)≤0,i∈{2, . . . ,m}.  (34)

To find the (unique) solution, we use Lagrange multipliers for convexoptimization: there exist λ_(i)≥0, λ∈{2, . . . , m}, such that

$\begin{matrix}{{{\frac{1}{2}{\nabla{\phi(y)}}} + {\lambda_{2}{\nabla{g_{2}(y)}}} + \ldots + {\lambda_{m}{\nabla{g_{m}(y)}}}} = 0} & (35) \\{{\forall{i \in {\left\{ {2,\ldots\mspace{11mu},m} \right\}\text{:}\mspace{14mu}\lambda_{i}{g_{i}(y)}}}} = 0} & (36) \\{\forall{i \in {{\left\{ {2,\ldots\mspace{11mu},m} \right\}\text{:}\mspace{14mu}{g_{i}(y)}} \leq 0.}}} & (37)\end{matrix}$

It follows from (35) that(y ₁−δ₁)+λ₂+ . . . +λ_(m)=0  (38)(y ₂−δ₂)−λ₂=0  (39). . .  (40)(y _(m)−δ_(m))−λ_(m)=0.  (41)

Then y₁≤δ₁ because λ₂, . . . , λ_(m)≥0. By substitution, one getsy ₁ +y ₂ + . . . +y _(m)=δ₁+δ₂+ . . . δ_(m).  (42)

Next, for each i∈{2, . . . , m}, one has y_(i)−δ_(i)=λ_(i)≥0,y_(i)−y₁≥0, and (36) implies,0=λ_(i) g _(i)(y)=(y _(i)−δ_(i))(y ₁ −y _(i)).  (43)This leads to either y_(i)=δ_(i)≥y₁ or y₁=y_(i)≥δ_(i). That meansy _(i)=max{y ₁,δ_(i)} for all i∈{2, . . . ,m}.  (44)

So (42) becomesδ₁+δ₂+ . . . +δ_(m) =y ₁+max{y ₁,δ₂}+ . . . +max{y ₁,δ_(m)}.  (45)Next, suppose k is the largest number in {1, . . . , m} such thatδ_(k)≤(δ₁+ . . . +δ_(k))/k.  (46)

The choice of k is always possible since at least (46) is true for k=1.Now setting δ_(m+1)=+∞, one claims that(δ₁+ . . . +δ_(k))/k<δ _(k+1).  (47)

In fact, if k=m, then (47) clearly holds. If k<m, then by the choice ofk, one has δ₁+ . . . +δ_(k)+δ_(k+1)<(k+1)δ_(k+1), which implies δ₁+ . .. +δ_(k)<kδ_(k+1). So (47) holds. Next, one can show thatδ_(k) ≤y ₁≤δ_(k+1).  (48)

Indeed, on the one hand, (46) and (45) imply

$\begin{matrix}{\mspace{20mu}{{{k\;\delta_{k}} + \delta_{k + 1} + \ldots + \delta_{m}} \leq {\left( {\delta_{1} + \ldots + \delta_{k}} \right) + \delta_{k + 1} + \ldots + \delta_{m}}}} & (49) \\{\mspace{20mu}{= {y_{1} + {\max\left\{ {y_{1},\delta_{2}} \right\}} + \ldots + {\max\left\{ {y_{1},\delta_{m}} \right\}}}}} & (50) \\{\mspace{20mu}{\leq {y_{1} + \underset{\underset{k - {1{terms}}}{︸}}{{\max\left\{ {y_{1},\delta_{k}} \right\}} + \ldots + {\max\left\{ {y_{1},\delta_{k}} \right\}}}}}} & (51) \\{\mspace{20mu}{{{+ \max}\left\{ {y_{1},\delta_{k + 1}} \right\}} + \ldots + {\max\left\{ {y_{1},\delta_{m}} \right\}}}} & (52) \\{= {y_{1} + {\left( {k - 1} \right)\max\left\{ {y_{1},\delta_{k}} \right)} + {\max\left\{ {y_{1},\delta_{k + 1}} \right\}} + \ldots + {\max{\left\{ {y_{1},\delta_{m}} \right\}.}}}} & (53)\end{matrix}$If y₁<δ_(k), then (4) implieskδ _(k) =y ₁+(k−1)max{y ₁,δ_(k) }<kδ _(k),  (54)which is a contradiction. Thus, y₁≥δ_(k). On the other hand, (47) and(45) implies

$\begin{matrix}{{{k\;\delta_{k + 1}} + \delta_{k + 1} + \ldots + \delta_{m}} \geq {\left( {\delta_{1} + \ldots + \delta_{k}} \right) + \delta_{k + 1} + \ldots + \delta_{m}}} & (55) \\{= {y_{1} + {\max\left\{ {y_{1},\delta_{2}} \right\}} + \ldots + {\max\left\{ {y_{1},\delta_{m}} \right\}}}} & (56) \\{= {y_{1} + \underset{\underset{k - {1{terms}}}{︸}}{y_{1} + \ldots + y_{1}} + {\max\left\{ {y_{1},\delta_{k + 1}} \right\}} + \ldots + {\max\left\{ {y_{1},\delta_{m}} \right\}}}} & (57) \\{\geq {{ky}_{1} + {\max\left\{ {y_{1},\delta_{k + 1}} \right\}} + \delta_{k + 2} + \ldots + {\delta_{m}.}}} & (58)\end{matrix}$If y₁>δ_(k+1), then (55)-(58) implies(k+1)δ_(k+1) ≥ky ₁+max{y ₁,δ_(k+1)}>(k+1)δ_(k+1),  (59)which is a contradiction. Thus, y₁≤δ_(k+1). Therefore, the claim (48) istrue. Then (45) impliesδ₁+δ₂+ . . . +δ_(m) =y ₁+max{y ₁,δ₂}+ . . . +max{y ₁,δ_(m) }=ky₁+δ_(k+1)+ . . . +δ_(m),   (60)which leads to y₁=(δ₁+ . . . +δ_(k))/k. Combining with (44) completesthe proof.

Edge Alignment Constraint

On the triangular mesh, the designer may want a constant slope on aparticular path, in which case one can say the path is “aligned”. Such apath is sometimes called a feature line. To formulate this constraint,suppose the feature line is given by adjacent points A₁, A₂, . . . ,A_(m) on the triangular mesh where A₁=(a_(i1), a_(i2), x_(i)). FIG. 13shows the birds-eye view of a feature line on a triangular mesh inaccordance with one or more embodiments of the invention. ForA_(i)A_(i+1), the length of its “shadow” on the xy-plane isδ_(i):=(a _(i+1,1) −a _(i,1))²+(a _(i+1,2) −a _(i,2))².  (61)Define alsot ₁:=0,t ₂:=δ₁ ,t ₃:=δ₁+δ₂ , . . . ,t _(m):=δ₁+ . . . +δ_(m−1).  (62)Then the alignment constraint is written as∀i∈{1, . . . ,m−2}:(x _(i+1) −x _(i))/(t _(i+1) −t _(i))=(x _(i+2) −x_(i+1))/(t _(i+2) −t _(i+1)).  (63)

Method Component 4 (Projection onto an Edge Alignment Constraint)

Suppose the points (a_(i1), a_(i2)), i∈{1, . . . , m} forms a path inR². Define (t₁, . . . , t_(m)) by (62) and define the convex setF:={(x ₁ , . . . ,x _(m))∈R ^(m)|((a _(i1) ,a _(i2) ,x_(i)))_(i∈{i, . . . ,m})satisfies(63)}⊂R ^(m).  (64)

Let z=(z₁, . . . , z_(m))∈R^(m). Then the projection P_(F)z is given by∀i ∈{1, . . . ,m}:(P _(F) Z)_(i)=ƒ(t _(i)).  (65)where ƒ(t):=α+βt is the least squares line for the points (t_(i),z_(i))_(i∈{1, . . . m}). More specifically, (α, β) is obtained from thelinear system

$\begin{matrix}{{\begin{bmatrix}m & {\sum\limits_{i = 1}^{m}t_{i}} \\{\sum\limits_{i = 1}^{m}t_{i}} & {\sum\limits_{i = 1}^{m}t_{i}^{2}}\end{bmatrix}\begin{bmatrix}\alpha \\\beta\end{bmatrix}} = \begin{bmatrix}{\sum\limits_{i = 1}^{m}z_{i}} \\{\sum\limits_{i = 1}^{m}{t_{i}z_{i}}}\end{bmatrix}} & (66)\end{matrix}$

To show this, first, F is convex since all constraints in (63) arelinear. Next, one considers the points (t_(i), z_(i)) in R². The problemis to find (x₁, . . . , x_(m)) such that the points (t_(i), x_(i)) arealigned and∥x−z∥ ² =E _(i=1) ^(m)(x _(i) −z _(i))² is minimized.  (67)

This is the least squares problem for the points (t_(i), z_(i)), shownin FIG. 14. So the least squares line is ƒ(t)=α+βt where (α, β)satisfies (66). Thus, the projection coordinates are (P_(F)(z))_(i)=α+βt_(i) for i∈{1, . . . , m}.

Surface Alignment Constraint

The designer may want to patch several adjacent triangles on the meshinto a single polygon, in which case one may say that these trianglesare “aligned”. This is equivalent to require all vertices of thetriangles to lie on the same plane. So the following may result.

Method Component 5 (Projection onto a Surface Alignment Constraint)

Let (a_(i1), a_(i2)), i∈{1, . . . , m} be a collection of points in R²that are not on the same line. DefineF:={(x ₁ , . . . ,x _(m))∈R ^(m)|{(a _(i1) ,a _(i2) ,x_(i))}_(i∈{1, . . . ,m})lie on the same plane}.   (69)Let z=(z₁, . . . , z_(m))∈R^(m). Then the projection P_(F)z is given by∀i∈{1, . . . ,m}:(P _(F) z)_(i)=ƒ(a _(i1) ,a _(i2)).  (70)where ƒ(t₁, t₂):=α+∈t₁+γt₂ is the least squares plane for the points(a_(i1), a_(i2), z_(i)), i∈{1, . . . , m}. More specifically, (α, β, γ)is the solution of the system

$\begin{matrix}{{\begin{bmatrix}m & \left\langle {e,a_{1}} \right\rangle & \left\langle {e,a_{1}} \right\rangle \\\left\langle {e,a_{1}} \right\rangle & {a_{1}}^{2} & {\left\langle {a_{1},a_{2}} \right\rangle} \\{\left\langle {e,a_{2}} \right\rangle} & \left\langle {a_{1},a_{2}} \right\rangle & {a_{2}}^{2}\end{bmatrix}\begin{bmatrix}\alpha \\\beta \\\gamma\end{bmatrix}} = {\begin{bmatrix}\left\langle {e,z} \right\rangle \\\left\langle {a_{1},z} \right\rangle \\\left\langle {a_{2},z} \right\rangle\end{bmatrix}.}} & (71)\end{matrix}$where e=(1, 1, . . . , 1), a₁=(a₁₁, a₂₁, . . . , a_(m1)), and a₂=(a₁₂a₂₂, . . . , a_(m2)).

To show this, first, F is a convex set. Let x=(x₁, . . . ,x_(m))=P_(F)Z, then x minimizes∥x−z∥ ²=Σ_(i=1) ^(m)(x _(i) −z _(i))²  (72)subject to the constraint that {(a_(i1), a_(i2),x_(i))}_(i∈{1, . . . , m}) lie on the same plane in R³. Thus, it isequivalent to find the least squares plane ƒ(t₁, t₂)=α+βt₁+γt₂ that fits{(a_(i1), a_(i2), z_(i))}_(i∈{1, . . . , m}). One defines e, a₁, a₂ asabove, then (α, β, γ) is the solution of (71). Since {(a_(i1),a_(i2))}_(i∈{1, . . . , m}) are not on the same line, (71) has a uniquesolution (α, β, γ). Finally, one computes x_(i)=α+βa_(i1)+γa_(i2) fori∈{1, . . . , m}.Projections onto General Surface Slope Constraints

Surface slope constraints are any requirement imposed on the slope of atriangle. In this section, a general setup is provided for projectionsonto such constraints.

Normal Vector and Surface Slope Constraint

Let (e₁, e₂, e₃) be the standard basis of R³. Given three pointsP₁=(P₁₁, P₁₂, h₁), P₂=(P₂₁, p₂₂, h₂), and P₃=(p₃₁, p₃₂, h₃) in R³, thenormal vector to the plane P₁P₂P₃ is the cross product

$\begin{matrix}{\overset{\rightarrow}{n} = {{\begin{bmatrix}{p_{11} - p_{31}} \\{p_{12} - p_{32}} \\{h_{1} - h_{3}}\end{bmatrix} \times \begin{bmatrix}{p_{21} - p_{31}} \\{p_{22} - p_{32}} \\{h_{2} - h_{3}}\end{bmatrix}} = {{{\text{:}\begin{bmatrix}a_{1} \\b_{1} \\t_{1}\end{bmatrix}} \times \begin{bmatrix}a_{2} \\b_{2} \\t_{2}\end{bmatrix}} = \begin{bmatrix}{{b_{1}t_{2}} - {b_{2}t_{1}}} \\{{a_{2}t_{1}} - {a_{1}t_{2}}} \\{{a_{1}b_{2}} - {a_{2}b_{1}}}\end{bmatrix}}}} & (73)\end{matrix}$The “shadows” of P₁, P₂, P₃ on xy-plane are (p₁₁, p₁₂), (p₂₁, p₂₂),(p₃₁, p₃₂), which are not on the same line. This implies a₁b₂−a₂b₁≠0. Soone can rescale and use

$\begin{matrix}{\overset{\rightarrow}{n} = {\left( {\frac{{b_{1}t_{2}} - {b_{2}t_{1}}}{{a_{1}b_{2}} - {a_{2}b_{1}}},{- \frac{{a_{1}t_{2}} - {a_{2}t_{1}}}{{a_{1}b_{2}} - {a_{2}b_{1}}}},1} \right).}} & (74)\end{matrix}$If one defines

${u = {{\frac{{b_{2}t_{1}} - {b_{1}t_{2}}}{{a_{1}b_{2}} - {a_{2}b_{1}}}\mspace{14mu}{and}\mspace{14mu} v} = \frac{{a_{1}t_{2}} - {a_{2}t_{1}}}{{a_{1}b_{2}} - {a_{2}b_{1}}}}},$thent ₁ =a ₁ u+b ₁ v,t ₂ =a ₂ u+b ₂ v, and {right arrow over(n)}=(−u,−v,1).  (75)Obviously, surface slope constraints depend solely on the normal vector{right arrow over (n)}. Thus, one can represent a general surface slopeconstraint in the formg(u,v)≤0.  (76)An important case of such constraints is the surface orientationconstraint below.

Definition 1 (surface orientation constraint) Let Δ be a triangle withnormal vector {right arrow over (n)}=(−u,−v,1)∈R³ as above. Let {rightarrow over (q)}=(q₁, q₂, q₃)∈R³ \{0} be a unit vector and θ be an anglein [0, π], the constraint

$\begin{matrix}{{{\angle\left( {\overset{\rightarrow}{n},\overset{\rightarrow}{q}} \right)} \leq \theta},{orequivalently},{{\cos\;{\angle\left( {\overset{\rightarrow}{n},\overset{\rightarrow}{q}} \right)}} = {\frac{{- {uq}_{1}} - {vq}_{2} + q_{3}}{\sqrt{u^{2} + v^{2} + 1}} \geq {\cos\;\theta}}},} & (77)\end{matrix}$is called the surface orientation constraint specified by the pair({right arrow over (q)}, θ). Below is a description of the generalframework for projection onto general surface slope constraints,followed by considerations of two special surface orientationconstraints: the surface maximum slope constraint and the surfaceoriented minimum slope constraint.

Projection onto a General Surface Slope Constraint

Consider a single triangle determined by three points Q₁=(p₁₁, p₁₂, w₁),Q₂=(p₂₁, p₂₂, w₂), and Q₃=(p₃₁, p₃₂, w₃) in R³. Without loss ofgenerality, one may assumew ₁ +w ₂ +w ₃=0.  (78)Projecting onto an orientation constraint is to find the new heights h₁,h₂, h₃ that minimizesϕ:=(h ₁ −w ₁)²+(h ₂ −w ₂)²+(h ₃ −w ₃)²  (79)such that the new points P₁=(p₁₁, p₁₂, h₁), P₂=(p₁₁, p₁₂, h₂), andP₃=(p₁₁, p₁₂, h₃) satisfy the given orientation constraint. Definet₁:=h₁−h₃ and t₂:=h₂−h₃, thenϕ=(t ₁ +h ₃ −w ₁)²+(t ₂ +h ₃ −w ₂)²+(h ₃ −w ₃)².  (80)Using the change of variables (75), one has:ϕ=(a ₁ u+b ₁ v+h ₃ −w ₁)²+(a ₂ u+b ₂ v+h ₃ −w ₂)²+(h ₃ −w ₃)².  (81)The new triangle P₁P₂P₃ needs to satisfy some given orientationconstraint g(u, v)≤0. So the projection reduces to

$\begin{matrix}{{\min\limits_{{({u,v,h_{3}})} \in R^{3}}\phi} = {\left( {{a_{1}u} + {b_{1}v} + h_{3} - w_{1}} \right)^{2} + \left( {{a_{2}u} + {b_{2}v} + h_{3} - w_{2}} \right)^{2} + \left( {h_{3} - w_{3}} \right)^{2}}} & (82) \\{\mspace{20mu}{{{subject}\mspace{14mu}{to}\mspace{14mu}{g\left( {u,v} \right)}} \leq 0.}} & (83)\end{matrix}$

Next, suppose that changing variables is necessary, for instance, (u, v)is replaced by the new variables (û, {circumflex over (v)}) under thelinear transformation

$\begin{matrix}{{\begin{bmatrix}\hat{u} \\\hat{v}\end{bmatrix}:={Q\begin{bmatrix}u \\v\end{bmatrix}}},{{where}\mspace{14mu} Q\mspace{14mu}{is}\mspace{14mu}{an}\mspace{14mu}{invertible}\mspace{14mu}{{matrix}.}}} & (84)\end{matrix}$

One may also define

$\begin{bmatrix}{\hat{a}}_{1} & {\hat{b}}_{1} \\{\hat{a}}_{2} & {\hat{b}}_{2}\end{bmatrix}:={\begin{bmatrix}a_{1} & b_{1} \\a_{2} & b_{2}\end{bmatrix}{Q^{- 1}.}}$Then

$\begin{matrix}{{{{{\hat{a}}_{1}\hat{u}} + {{\hat{b}}_{1}\hat{v}}} = {{\begin{bmatrix}{\hat{a}}_{1} & {\hat{b}}_{1}\end{bmatrix}\begin{bmatrix}\hat{u} \\\hat{v}\end{bmatrix}} = {{\begin{bmatrix}a_{1} & b_{1}\end{bmatrix}Q^{- 1}{Q\begin{bmatrix}u \\v\end{bmatrix}}} = {{a_{1}u} + {b_{1}v}}}}},{and}} & (85) \\{\mspace{20mu}{{{{\hat{a}}_{2}\hat{u}} + {{\hat{b}}_{2}\hat{v}}} = {{a_{2}u} + {b_{2}{v.}}}}} & (86)\end{matrix}$So (82)-(83) becomes

$\begin{matrix}{{\min\limits_{{({\hat{u},\hat{v},h_{3}})} \in R^{3}}\phi} = {\left( {{{\hat{a}}_{1}\hat{u}} + {{\hat{b}}_{1}\hat{v}} + h_{3} - w_{1}} \right)^{2} + \left( {{{\hat{a}}_{2}\hat{u}} + {{\hat{b}}_{2}\hat{v}} + h_{3} - w_{2}} \right)^{2} + \left( {h_{3} - w_{3}} \right)^{2}}} & (87) \\{\mspace{20mu}{{{{subject}\mspace{14mu}{to}\mspace{14mu}{\hat{g}\left( {\hat{u},\hat{v}} \right)}}:={{g\left( {u,v} \right)} \leq 0}},}} & (88)\end{matrix}$Therefore, one can treat (87)-(88) as (82)-(83) with new respectivecoefficients â_(i), {circumflex over (b)}_(i)

Next, one can simplify the model (82)-(83). Since (83) does not involveh₃, one can convert problem (82)-(83) into two variables (u, v) asfollows: first, set the derivative of ϕ with respect to h₃ to zero∇_(h) ₃ ϕ(u,v,h ₃)=2(a ₁ u+b ₁ v+h ₃ −w ₁)+2(a ₂ u+b ₂ v+h ₃ −w ₂)+2(h ₃−w ₃)=0.(89)Then one obtains

$\begin{matrix}{h_{3} = {- {{\frac{1}{3}\left\lbrack {{\left( {a_{1} + a_{2}} \right)u} + {\left( {b_{1} + b_{2}} \right)v}} \right\rbrack}.}}} & (90)\end{matrix}$Substituting into (82) and multiplying by 9, one obtains

$\begin{matrix}{{9\phi} = \left\lbrack {{\left( {{2a_{1}} - a_{2}} \right)u} + {\left( {{2b_{1}} - b_{2}} \right)v} - {3w_{1}}} \right\rbrack^{2}} & (91) \\{+ \left\lbrack {{\left( {{2a_{2}} - a_{1}} \right)u} + {\left( {{2b_{2}} - b_{1}} \right)v} - {3w_{2}}} \right\rbrack^{2}} & (92) \\{+ \left\lbrack {{\left( {a_{1} + a_{2}} \right)u} + {\left( {b_{1} + b_{2}} \right)v} - {3\left( {w_{1} + w_{2}} \right)}} \right\rbrack^{2}} & (93) \\{= {{\text{:}\left( {{6a_{1}^{2}} + {6a_{2}^{2}} - {6a_{1}a_{2}}} \right)u^{2}} + {\left( {{6b_{1}^{2}} + {6b_{2}^{2}} - {6b_{1}b_{2}}} \right)v^{2}}}} & (94) \\{{+ 6}\left( {{2a_{1}b_{1}} + {2a_{2}b_{2}} - {a_{1}b_{2}} - {a_{2}b_{1}}} \right){uv}} & (95) \\{{{- 18}\left( {{a_{1}w_{1}} + {a_{2}w_{2}}} \right)u} - {18\left( {{b_{1}w_{1}} + {b_{2}w_{2}}} \right)v} + {Z.}} & (96)\end{matrix}$where Z is some constant. So, by definingA=2(a ₁ ² +a ₂ ² −a ₁ a ₂)>0,  (97)B=2(b ₁ ² +b ₂ ² −b ₁ b ₂)>0,  (98)C=2a ₁ b ₁+2a ₂ b ₂ −a ₁ b ₂ −a ₂ b ₁,  (99)w _(a)=3(a ₁ w ₁ +a ₂ w ₂),  (100)w _(b)=3(b ₁ w ₁ +b ₂ w ₂),  (101)one has 9ϕ=3Au²+3Bv²+6Cuv−6w_(a)u−6w_(b)v+Z. Thus, (82)-(83) reduces to

$\begin{matrix}{{\min\limits_{{({u,v})} \in R^{2}}{\phi\left( {u,v} \right)}} = {{\frac{A}{2}u^{2}} + {Cuv} + {\frac{B}{2}v^{2}} - {w_{a}u} - {w_{b}v}}} & (102) \\{{{subjectto}\mspace{14mu}{g\left( {u,v} \right)}} \leq 0.} & (103)\end{matrix}$As long as one can find the solution (u, v), one can obtain theprojection (h₁, h₂, h₃) from (90), (75), and (73). More specifically,

$\begin{matrix}{{h_{3} = {- {\frac{1}{3}\left\lbrack {{\left( {a_{1} + a_{2}} \right)u} + {\left( {b_{1} + b_{2}} \right)v}} \right\rbrack}}},} & (104) \\{{h_{1} = {{a_{1}u} + {b_{1}v} + h_{3}}},} & (105) \\{h_{2} = {{a_{2}u} + {b_{2}v} + {h_{3}.}}} & (106)\end{matrix}$Due to (85)-(86), if new variables (û, {circumflex over (v)}) are used,then one also has (104)-(106) in which (a_(i), b_(i), u, v) is replacedby (â_(i), {circumflex over (b)}_(i), û, {circumflex over (v)}).Projection onto Surface Maximum Slope Constraint

Let P₁P₂P₃ be a triangle in R³ with normal vector {right arrow over(n)}=(−u,−v,1) obtained as in the normal vector and surface slopeconstraint section above. In certain cases, it is required that theangle between {right arrow over (n)} and a given direction {right arrowover (q)} must not exceed a given threshold. For example, suppose P₁P₂P₃represents the desired ground, that cannot be too steep with respect togravity, i.e., the slope of the plane P₁P₂P₃ must not exceed a thresholds:=s_(max)∈R₊. Then the angle between {circumflex over (n)} ande₃=(0,0,1) must satisfy ∠({circumflex over (n)},e₃)≤tan⁻¹ (s), which isshown as cone 1502 in FIG. 15. Accordingly, FIG. 15 illustrates amaximum angle for a normal vector of a ground plane with respect togravity in accordance with one or more embodiments of the invention.Thus, whenever the normal vector {right arrow over (n)} is outside thecone 1502, the constraint is violated. The inequality ∠({right arrowover (n)}, e₃)≤tan⁻¹ (s) is equivalent tocos ∠({right arrow over (n)},e ₃)=

{right arrow over (n)},e ₃

/∥{right arrow over (n)}∥≥cos(tan⁻¹(s))=1/√{square root over (1+s²)}.  (107)

Definition 2 (surface maximum slope constraint) One calls (107) thesurface maximum slope constraint with maximum slope s. This is a specialcase of surface orientation constraint where ({right arrow over (q)},θ)=(e₃, tan⁻¹ (s)). From (107), one may conclude

$\begin{matrix}\left. {\frac{1}{\sqrt{u^{2} + v^{2} + 1}} \geq \frac{1}{\sqrt{1 + s^{2}}}}\Leftrightarrow{{u^{2} + v^{2} - s^{2}} \leq 0.} \right. & (108)\end{matrix}$Using the general model (102)-(103) with the constants defined by(97)-(101), projection onto the surface maximum slope constraint reducesto

$\begin{matrix}{{\min\limits_{{({u,v})} \in R^{2}}\phi} = {{\frac{A}{2}u^{2}} + {Cuv} + {\frac{B}{2}v^{2}} - {w_{a}u} - {w_{b}v}}} & (109) \\{{{{{subject}\mspace{14mu}{to}\mspace{14mu} u^{2}} + v^{2} - s^{2}} \leq 0},} & (110)\end{matrix}$Rescaling by

$\left. \left( {u,v,w_{a},w_{b}} \right)\leftarrow\left( {\frac{u}{s},\frac{v}{s},\frac{w_{a}}{s},\frac{w_{b}}{s}} \right) \right.,$one can obtain an equivalent problem

$\begin{matrix}{{\min\limits_{{({u,v})} \in R^{2}}{\phi\left( {u,v} \right)}} = {{\frac{A}{2}u^{2}} + {Cuv} + {\frac{B}{2}v^{2}} - {w_{a}u} - {w_{b}u}}} & (111) \\{{{subject}\mspace{14mu}{to}\mspace{14mu}{g\left( {u,v} \right)}}:={{u^{2} + v^{2} - 1} \leq 0.}} & (112)\end{matrix}$

Hence, the projection to the surface maximum slope constraint will findnew elevations for the triangle vertices, such that the angle betweenthe z-axis and the normal vector of the triangle surface does not exceedtan⁻¹ (s) (i.e. the normal vector lays in the cone 1502 depicted in FIG.15), and that equation (79) is minimized.

This problem may be solved directly by several ways including numericalmethods. One exemplary method is presented herein, in which thewell-known Lagrange multipliers, also known as Karush-Kuhn-Tucker (KKT)conditions [19, 20], and Ferrari's method for quartic equations [16] areused. First, one observes the following:

1. ϕ(u, v) is a strictly convex quadratic surface.

2. For each value η≥0, the level set ϕ(u, v)=η is an ellipse in R² whosecenter is the vertex (u₀, v₀) of ϕ(u, v), which satisfiesAu ₀ +Cv ₀ =w _(a)  (113)Cu ₀ +Bv ₀ =w _(b).  (114)

3. This is minimizing a strictly convex function over a closed, bounded,convex set. Thus, there exists a unique solution.

One may consider two cases:

Case 1: g(u₀, v₀)≤0. Then (u₀, v₀) is the solution of (111)-(112).

Case 2: g(u₀, v₀)>0. Then the solution of (111)-(112) is the tangentpoint of the circle g(u, v)=0 and the smallest elliptical level set ϕ(u,v)=r that intersects the circle (see FIG. 16). Accordingly, FIG. 16illustrates an optimal solution based on tangent points between a circleand level sets. To find the maximizers/minimizers, one considers thefollowing Lagrange multipliers system

$\begin{matrix}{{{{\nabla_{u}{\phi\left( {u,v} \right)}} + {\frac{\lambda}{2}{\nabla_{u}{g\left( {u,v} \right)}}}} = 0},} & (115) \\{{{{\nabla_{u}{\phi\left( {u,v} \right)}} + {\frac{\lambda}{2}{\nabla_{v}{g\left( {u,v} \right)}}}} = 0},} & (116) \\{{u^{2} + v^{2} - 1} = 0.} & (117)\end{matrix}$

Using the geometry of R² (see FIG. 16) and the theory of Lagrangemultipliers, one can obtain the below properties.

Consider problem (111)-(112) and its Lagrange multipliers system(115)-(117). Suppose the solution (u₀, v₀) of (113)-(114) satisfiesg(u₀, v₀)>0. Then

1. (115)-(117) yields a unique solution (u₁, v₁, λ₁) such that λ₁>0; andin this case, (u₁, v₁) is the unique solution of (111)-(112).

2. (115)-(117) yields at least one solution (u₂, v₂, λ₂) such that λ₂<0.

Therefore, one seeks a solution (u, v, λ) of (115)-(117) with λ>0.First, (115) and (116) yield(A+λ)u+Cv=w _(a)  (118)Cu+(B+λ)v=w _(b).  (119)

Using the Cauchy-Schwarz inequality, one hasAB=[a ₁ ² +a ₂ ²+(a ₁ −a ₂)²][b ₁ ² +b ₂ ²+(b ₁ −b ₂)²](120)≥[a ₁ b ₁ +a ₂ b ₂+(a ₁ −a ₂)(b ₁ −b ₂)]² =C ².  (121)

Since A, B, λ>0, the determinant of (118)-(119) is(A+λ)(B+λ)−C ² >AB−C ²>0.  (122)This implies the system (118)-(119) has a unique solution

$\begin{matrix}{{u = {\frac{{w_{a}\left( {B + \lambda} \right)} - {w_{b}C}}{{\left( {A + \lambda} \right)\left( {B + \lambda} \right)} - C^{2}} = \frac{{\lambda\; w_{a}} + {w_{a}B} - {w_{b}C}}{\lambda^{2} + {\left( {A + B} \right)\lambda} + {AB} - C^{2}}}},} & (123) \\{{v = {\frac{{w_{b}\left( {A + \lambda} \right)} - {w_{a}C}}{{\left( {A + \lambda} \right)\left( {B + \lambda} \right)} - C^{2}} = \frac{{\lambda\; w_{b}} + {w_{b}A} - {w_{a}C}}{\lambda^{2} + {\left( {A + B} \right)\lambda} + {AB} - C^{2}}}},} & (124)\end{matrix}$

Next, substituting u and v into (117) yields

$\begin{matrix}\left\lbrack {\lambda^{2} + {\left( {A + B} \right)\lambda} + {AB} - C^{2}} \right\rbrack^{2} & (125) \\{= {\left\lbrack {{w_{a}\lambda} + \left( {{w_{a}B} - {w_{b}C}} \right)} \right\rbrack^{2} + \left\lbrack {{w_{b}\lambda} + \left( {{w_{b}A} - {w_{a}C}} \right)} \right\rbrack^{2}}} & (126) \\{= {{\left( {w_{a}^{2} + w_{b}^{2}} \right)\lambda^{2}} + {2\left( {{w_{a}^{2}B} + {w_{b}^{2}A} - {2w_{a}w_{b}C}} \right)\lambda}}} & (127) \\{{+ \left( {{w_{a}B} - {w_{b}C}} \right)^{2}} + {\left( {{w_{b}A} - {w_{a}C}} \right)^{2}.}} & (128)\end{matrix}$

Defining the constants accordingly, one may rewrite as(λ² +R ₁ λ+R ₂)² =R ₃λ²+2R ₄ λ+R ₅.  (129)

One can simply solve this equation by the classic Ferrari's method forquartic equations [9]. However, since there is a unique positive λ, onewill find its explicit formula following Ferrari's technique. A realvariable y is introduced(λ² +R ₁ λ+R ₂ +y)² =R ₃λ²+2R ₄ λ+R ₅+2(λ² +R ₁ λ+R ₂)y+y ²  (130)=(R ₃+2y)λ²+2(R ₄ +R ₁ y)+R ₅+2R ₂ y+y ².  (131)

One can choose y such that the right hand side is a perfect square in λ.Thus, the right hand side must have zero discriminant, i.e.,[R ₄ +R ₁ y]²−(R ₃+2y)(R ₅+2R ₂ y+y ²)=0

  (132)−2y ³+[R ₁ ² −R ₃−4R ₂]y ²+[2R ₁ R ₄−2R ₂ R ₃−2R ₅]y+R ₄ ² −R ₃ R ₅=0.  (133)

This is a cubic equation in y, so Cardano's method [9] can be used tofind one real solution y₀. Then (130)-(131) becomes

$\begin{matrix}{\left( {\lambda^{2} + {R_{1}\lambda} + R_{2} + y_{0}} \right)^{2} = {\left( {R_{3} + {2y_{0}}} \right){\left( {\lambda + \frac{R_{4} + {R_{1}y_{0}}}{R_{3} + {2y_{0}}}} \right)^{2}.}}} & (134)\end{matrix}$From the discussion above, this equation must have at least one positiveand one negative solutions. Next, one solves for the (unique) positiveΔ:

Case 2a: R₃+2y₀<0. Then

$\lambda = {- \frac{R_{4} + {R_{1}y_{0}}}{R_{3} + {2y_{0}}}}$is the unique value of λ, which is a contradiction. Thus, this casecannot happen.

Case 2b. R₃+²y₀=0. Then

$y_{0} = {- \frac{R_{3}}{2}}$andλ² +R ₁ λ+R ₂ +y ₀=0.  (135)So one has

$\begin{matrix}{\lambda = {{\frac{1}{2}\left( {{- R_{1}} \pm \sqrt{R_{1}^{2} - {4\left( {R_{2} - \frac{R_{3}}{2}} \right)}}} \right)} = {\frac{1}{2}{\left( {{- R_{1}} \pm \sqrt{R_{1}^{2} - {4R_{2}} + {2R_{3}}}} \right).}}}} & (136)\end{matrix}$Since there is only one positive λ, one takes

$\begin{matrix}{\lambda = {\frac{1}{2}{\left( {{- R_{1}} + \sqrt{R_{1}^{2} - {4R_{2}} + {2R_{3}}}} \right).}}} & (137)\end{matrix}$

Case 2c: R₃+2y₀>0. Define r:=√{square root over (R₃+2y₀)}. From (134),

$\begin{matrix}{\left( {\lambda^{2} + {R_{1}\lambda} + R_{2} + y_{0}} \right)^{2} = \left( {{\lambda\sqrt{R_{3} + {2y_{0}}}} + \frac{R_{4} + {R_{1}y_{0}}}{\sqrt{R_{3} + {2y_{0}}}}} \right)^{2}} & (138) \\{= {\left( {{\lambda\; r} + \frac{R_{4} + {R_{1}y_{0}}}{r}} \right)^{2}.}} & (139)\end{matrix}$This leads to two quadratic equations in λ

$\begin{matrix}{{\lambda^{2} + {\left( {R_{1} \pm r} \right)\lambda} + \left( {R_{2} + {y_{0} \pm \frac{R_{4} + {R_{1}y_{0}}}{r}}} \right)} = 0.} & (140)\end{matrix}$

Since there is only one positive λ, one must have R₄+R₁y₀≠0. Otherwise,(140) consists of two quadratic equations that have constant term R₂+y₀;thus, will yield an even number (possibly none) of positive solutions,which is a contradiction. Moreover, one must take the equation withnegative constant term. So one sets r←r·sgn(R₄+R₁y₀) and takes only theequation

$\begin{matrix}{{\lambda^{2} + {\left( {R_{1} - r} \right)\lambda} + \left( {R_{2} + y_{0} - \frac{R_{4} + {R_{1}y_{0}}}{r}} \right)} = 0.} & \left. (141) \right)\end{matrix}$It follows that the positive Δ is

$\begin{matrix}{\lambda = {{\frac{1}{2}\left\lbrack {r - R_{1} + \sqrt{\left( {r - R_{1}} \right)^{2} - {4\left( {R_{2} + y_{0} - \frac{R_{4} + {R_{1}y_{0}}}{r}} \right)}}} \right\rbrack}.}} & (142)\end{matrix}$

Next, one obtains u and v from (123)-(124) and then rescales variables(u, v)←(su, sv). Finally, one can use (104)-(106) to obtain (h₁, h₂,h₃).

Projection onto Surface Oriented Minimum Slope Constraint

Nonconvex Surface Minimum Constraint

Let P₁P₂P₃ be a triangle with the normal vector {right arrow over(n)}=(−u, −v, 1) as described above. In some cases, it is required thatthe angle ∠({right arrow over (n)}, {right arrow over (q)})≥α for somegiven vector {right arrow over (q)} and number α.

This happens, for example, if P₁P₂P₃ must have a slope at leasts:=s_(min) ∈R₊. In this case, the angle between {right arrow over (n)}and e₃=(0,0,1) satisfies∠({right arrow over (n)},e ₃)≥tan⁻¹(s).  (143)

If one considers FIG. 15 again, in this case, the normal vector needs tobe outside of the cone 1502. The inequality ∠({right arrow over (n)},e₃)≥tan⁻¹ (s) is equivalent to

$\begin{matrix}{{\cos\;{\angle\left( {\overset{->}{n},e_{3}} \right)}} = {{\frac{\left\langle {\overset{->}{n},e_{3}} \right\rangle}{\overset{->}{n}} \leq {\cos\;\left( {\tan^{- 1}(s)} \right)}} = {\frac{1}{\sqrt{1 + s^{2}}}.}}} & (144)\end{matrix}$One can refer to (143) or (144) as the surface minimum slope constraintwith minimum slope s. From (144), one has

$\begin{matrix}\left. {\frac{1}{\sqrt{u^{2} + v^{2} + 1}} \leq \frac{1}{\sqrt{1 + s^{2}}}}\Leftrightarrow{{u^{2} + v^{2} - s^{2}} \geq 0.} \right. & (145)\end{matrix}$This is clearly a nonconvex constraint in (u, v). Despite the projectiononto this constraint that is still available, nonconvexity may preventiterative methods from convergence. Therefore, one will find a convexreplacement for this constraint which can still guarantee (143)-(144).

The Surface Oriented Minimum Slope Constraint

(143) implies the angle between {right arrow over (n)} and the xy-planesatisfies∠({right arrow over (n)},xy−plane)≤(π/2)−tan⁻¹(s).  (146)In some cases, it is not unreasonable to align the plane P₁P₂P₃ toward agiven location/direction. For instance, in civil engineering drainage,the designer may want to direct the water to certain drains or inlets.Suppose one wants P₁P₂P₃ to incline toward a unit direction {right arrowover (q)}=(q₁, q₂, 0)∈R³ (see FIG. 17). In this regard, FIG. 17illustrates a desired orientation ({right arrow over (q)}) of a normalvector (n) for a triangle plane. Then one can fulfill (146) by requiring

${{\angle\left( {\overset{->}{n},\overset{->}{q}} \right)} \leq {\frac{\pi}{2} - {\tan^{- 1}(s)}}},$i.e.,

$\begin{matrix}{{\cos\;{\angle\left( {\overset{->}{n},\overset{->}{q}} \right)}} = {{\frac{\left\langle {\overset{->}{n},\overset{->}{q}} \right\rangle}{\overset{->}{n}} \geq {\cos\;\left( {\frac{\pi}{2} - {\tan^{- 1}(s)}} \right)}} = \frac{s}{\sqrt{1 + s^{2}}}}} & (147)\end{matrix}$and arrive at the following definition.

Definition 3 [surface oriented minimum slope constraint]. One calls(147) the surface oriented minimum slope constraint specified by ({rightarrow over (q)}, s), the unit vector {right arrow over (q)} and theminimum slope s. This is a special case of the surface orientationconstraint where

$\begin{matrix}{\left( {\overset{->}{q},\theta} \right) = {\left( {\overset{->}{q},{\frac{\pi}{2} - {\tan^{- 1}(s)}}} \right).}} & \;\end{matrix}$Substituting {right arrow over (n)}=(−u, −v, 1) and {right arrow over(q)}=(q₁, q₂, 0) into (147), one has

$\begin{matrix}{{\frac{{{- q_{1}}u} - {q_{2}v}}{\sqrt{u^{2} + v^{2} + 1}} \geq \frac{s}{\sqrt{1 + s^{2}}}},} & (148)\end{matrix}$which is equivalent to

$\begin{matrix}{{{{q_{1}u} + {q_{2}v}} < 0},} & (149) \\{{u^{2} + v^{2} + 1 - {\frac{1 + s^{2}}{s^{2}}\left( {{q_{1}u} + {q_{2}v}} \right)^{2}}} \leq 0.} & (150)\end{matrix}$This constraint not only guarantees surface minimum slope, but alsogenerates a convex set in (u, v). A visualization of this finding isgiven in FIG. 18. In other words, FIG. 18 illustrates a visualization ofa surface minimum slope constraint in accordance with one or moreembodiments of the invention. The constraint requires that the normalvector lays on or inside the half cone 1802, whereas the directionvector q is the axis of the full cone.

Projection onto Surface Oriented Minimum Slope Constraint

By the model (82)-(83), the projection onto the surface oriented minimumslope constraint is

$\begin{matrix}{{\min\limits_{{({u,v})} \in R^{2}}{\phi\left( {u,v,h_{3}} \right)}} = {\left( {{a_{1}u} + {b_{1}v} + h_{3} - w_{1}} \right)^{2} + \left( {{a_{2}u} + {b_{2}v} + h_{3} - w_{2}} \right)^{2} + \left( {h_{3} - w_{3}} \right)^{2}}} & (151) \\{\mspace{79mu}{{{subject}\mspace{14mu}{to}\mspace{14mu}(149)} - {(150).}}} & (152)\end{matrix}$

Again, one solves (151)-(152) directly using Lagrange multipliers andFerrari's method for quartic equations. First, define

${Q:=\begin{bmatrix}q_{1} & q_{2} \\q_{2} & {- q_{1}}\end{bmatrix}},$then Q²=Id since∥{right arrow over (q)}∥=q₁ ²+q₂ ²=1. So Q=Q^(T)=Q⁻¹.Next, the variables are changed

$\begin{matrix}{\begin{bmatrix}\hat{u} \\\hat{v}\end{bmatrix}:={{Q\begin{bmatrix}u \\v\end{bmatrix}} = {\left. \begin{bmatrix}{q_{1}u} & {q_{2}v} \\{q_{2}u} & {{- q_{1}}v}\end{bmatrix}\Leftrightarrow\begin{bmatrix}u \\v\end{bmatrix} \right. = {{Q\begin{bmatrix}\hat{u} \\\hat{v}\end{bmatrix}} = {\begin{bmatrix}{q_{1}\hat{u}} & {q_{2}\hat{v}} \\{q_{2}\hat{u}} & {{- q_{1}}\hat{v}}\end{bmatrix}.}}}}} & (153)\end{matrix}$Then (150) becomes

$\begin{matrix}{{\left( {{q_{1}\hat{u}} + {q_{2}\hat{v}}} \right)^{2} + \left( {{q_{2}\hat{u}} - {q_{1}\hat{v}}} \right)^{2} + 1 - {\left( {1 + \frac{1}{s^{2}}} \right){\hat{u}}^{2}}} = {{{\left( {- \frac{1}{s^{2}}} \right){\hat{u}}^{2}} + {\hat{v}}^{2} + 1} \leq 0.}} & (154)\end{matrix}$Thus, the constraint (149)-(150) becomes

$\begin{matrix}{{\hat{u} < 0},} & (155) \\{{{{\hat{v}}^{2} - \frac{{\hat{u}}^{2}}{s^{2}} + 1} \leq 0},} & (156)\end{matrix}$

Following the projections onto a general surface slope constraintdescribed above, the coefficients are changed (without relabeling)

$\begin{matrix}\left. \begin{bmatrix}a_{1} & a_{2} \\b_{1} & b_{2}\end{bmatrix}\leftarrow{Q\begin{bmatrix}a_{1} & a_{2} \\b_{1} & b_{2}\end{bmatrix}} \right. & (157)\end{matrix}$and the following problem is obtained

$\begin{matrix}{{\min\limits_{{({\hat{u},\hat{v}})} \in R^{2}}{\phi\left( {\hat{u},\hat{v}} \right)}} = {{\frac{A}{2}{\hat{u}}^{2}} + {C\hat{u}\hat{v}} + {\frac{B}{2}{\hat{v}}^{2}} - {w_{a}\hat{u}} - {w_{b}\hat{v}}}} & (158) \\{{subject}\mspace{14mu}{to}\mspace{14mu}\left\{ {\begin{matrix}{{\hat{u} < 0},} \\{{{\hat{v}}^{2} - \frac{{\hat{u}}^{2}}{s^{2}} + 1} \leq 0}\end{matrix}.} \right.} & (159)\end{matrix}$where A, B, C, w_(a), w_(b) are defined below using the new a₁, a₂, b₁,b₂ from (157)A=2(a ₁ ² a ₂ ² −a ₁ a ₂)>0,  (160)B=2(b ₁ ² +b ₂ ² −b ₁ b ₂)>0,  (161)C=2a ₁ b ₁+2a ₂ b ₂ −a ₁ b ₂ −a ₂ b ₁,  (162)w _(a)=3(a ₁ w ₁ +a ₂ w ₂),  (163)w _(b)=3(b ₁ w ₁ +b ₂ w ₂).  (164)Changing variables by

$\left. \left( {u,v,A,B,w_{b}} \right)\leftarrow\left( {\frac{\hat{u}}{s},\hat{v},{sA},\frac{B}{s},\frac{w_{b}}{s}} \right) \right.,$then (158)-(159) becomes

$\begin{matrix}{{\min\limits_{{({\hat{u},\hat{v}})} \in R^{2}}{\phi\left( {u,v} \right)}} = {{\frac{A}{2}u^{2}} + {Cuv} + {\frac{B}{2}v^{2}} - {w_{a}u} - {w_{b}v}}} & (165) \\{{subjectto}\mspace{14mu}\left\{ {\begin{matrix}{{{g_{1}\left( {u,v} \right)}:={u < 0}},} \\{{g_{2}\left( {u,v} \right)}:={{v^{2} - u^{2} + 1} \leq 0}}\end{matrix}.} \right.} & (166)\end{matrix}$

Hence, the projection to the surface oriented minimum slope constraintwill find new elevations for the triangle vertices, such that the anglebetween the normal vector of the triangle and the z-axis is at leasttan⁻¹(s) and the angle between the normal vector and the directionalvector q is at most

${\frac{\pi}{2} - {\tan^{- 1}(s)}},$and that equation (151) is minimized.

Before solving the above problem, one can observe that

1. The constraint (166) is one side of a hyperbola, thus, convex.

2. The objective ϕ(u, v) in (165) is a nonnegative strictly convexquadratic surface.

3. For each value η≥0, the level set ϕ(u, v)=η is an ellipse in R² whosecenter is the vertex (u₀, v₀) of ϕ(u, v), which satisfiesAu ₀ +Cv ₀ =w _(a)  (167)Cu ₀ +Bv ₀ =w _(b).  (168)

Therefore, (165)-(166) has a unique solution. To solve it, two cases areconsidered:

Case 1: (u₀, v₀) is feasible, i.e., it lies inside or on the left branch{(u, v)∈R²|g₂(u, v)=0, u<}. Then it is the unique solution of(165)-(166) because the optimal objective value is zero. FIG. 19illustrates a feasible set region for a normal vector with respect to aSurface Oriented Minimum Slope Constraint in accordance with one or moreembodiments of the invention.

Case 2: (u₀, v₀) is not feasible. Let (u, v) be a solution to(165)-(166). By Lagrange multipliers, there exists λ>0 such that

$\begin{matrix}{{{{\nabla_{u}{\phi\left( {u,v} \right)}} + {\frac{\lambda}{2}{\nabla_{u}{g_{2}\left( {u,v} \right)}}}} = 0},} & (169) \\{{{{\nabla_{v}{\phi\left( {u,v} \right)}} + {\frac{\lambda}{2}{\nabla_{v}{g_{2}\left( {u,v} \right)}}}} = 0},} & (170) \\{{{g_{2}\left( {u,v} \right)} = {{v^{2} - u^{2} + 1} = 0}},} & (171) \\{{g_{1}\left( {u,v} \right)} = {u < 0.}} & (172)\end{matrix}$

Consider the system of three equations (169)-(171). Using the geometryof R² (see FIGS. 20 and 21), one may conclude that it has at least twosolutions. To describe its solutions, let (u₀, v₀) be defined by(167)-(168), then

1. If g₂(u₀, v₀)>0, then (169)-(171) has a unique solution (u₁, v₁, λ₁)with λ₁>0 and u₁<0, and a unique solution (u₂, v₂, λ₂) with λ₂>0 andu₂>0.

2. If g₂(u₀, v₀)≤0 and g₁(u₀, v₀)=u>0, then (169)-(171) has a uniquesolution (u₁, v₁, λ₁) with λ₁>0. Thus, u₁<0. In addition, (169)-(171)also has at least one solution (u₂, v₂, λ₂) with λ₂<0.

3. In both cases 1 and 2, (u₁, v₁, λ₁) is the unique solution of thewhole Lagrange system (169)-(172), thus, (u₁, v₁) is the unique solutionof (165)-(166).

To prove the above conclusion, suppose (u, v, λ) is a solution of(169)-(171), then (u, v) is the tangent point of the hyperbola v²−u²+1=0and an elliptical level set of ϕ(u, v).

FIG. 20 shows that there are exactly two tangent points (u, v) withpositive Lagrange multipliers λ, one of which has u<0 and the other hasu>0 in accordance with one or more embodiments of the invention. Thisjustifies conclusion 1.

FIG. 21 shows that there is exactly one tangent point (u, v) with apositive Lagrange multiplier λ and u>0, and at least one tangent pointwith a nonpositive Lagrange multiplier in accordance with one or moreembodiments of the invention. This justifies conclusion λ.

Now if (u, v) is a solution to (165)-(166), then there exists λ>0 suchthat (u, v, λ) is a solution to (169)-(171). So, conclusion 3 followsfrom conclusions 1 and 2.

Next, one can show how to find the solution in Case 2. First, write(169)-(170) as(A−λ)u+Cv=w _(a),  (173)Cu+(B+λ)v=w _(b).  (174)Suppose for now that the determinant (A−λ)(B+λ)−C²≠0. Then

$\begin{matrix}{{{{u = \frac{{w_{a}\left( {B + \lambda} \right)} - {w_{b}C}}{{\left( {A - \lambda} \right)\left( {B + \lambda} \right)} - C^{2}}}\quad} = \frac{\left( {{\lambda\; w_{a}} + {w_{a}B} - {w_{b}C}} \right)}{{- \lambda^{2}} + {\left( {A - B} \right)\lambda} + \left( {{AB} - C^{2}} \right)}},} & (175) \\{{{{v = \frac{{w_{b}\left( {A - \lambda} \right)} - {w_{a}C}}{{\left( {A - \lambda} \right)\left( {B + \lambda} \right)} - C^{2}}}\quad} = \frac{{{- \lambda}\; w_{b}} + \left( {{w_{b}A} - {w_{a}C}} \right)}{{- \lambda^{2}} + {\left( {A - B} \right)\lambda} + \left( {{AB} - C^{2}} \right)}},} & (176)\end{matrix}$Substituting into (171), one can derive

$\begin{matrix}{1 = {\left( \frac{{\lambda\; w_{a}} + {w_{a}B} - {w_{b}C}}{\lambda^{2} + {\left( {B - A} \right)\lambda} + \left( {C^{2} - {AB}} \right)} \right)^{2} - {\left( \frac{{\lambda\; w_{b}} + \left( {{w_{a}C} - {w_{b}A}} \right)}{\lambda^{2} + {\left( {B - A} \right)\lambda} + \left( {C^{2} - {AB}} \right)} \right)^{2}.}}} & (177)\end{matrix}$This implies

$\begin{matrix}\begin{matrix}\left\lbrack {\lambda^{2} + {\left( {B - A} \right)\lambda} + \left( {C^{2} - {AB}} \right)} \right\rbrack^{2} \\{= {\left( {{\lambda\; w_{a}} + {w_{a}B} - {w_{b}C}} \right)^{2} - \left\lbrack {{\lambda\; w_{b}} + \left( {{w_{a}C} - {w_{b}A}} \right)} \right\rbrack^{2}}} \\{= {{\left( {w_{a}^{2} - w_{b}^{2}} \right)\lambda^{2}} + {2\left( {{w_{a}^{2}B} + {w_{b}^{2}A} - {2w_{a}w_{b}C}} \right)\lambda}}} \\{{+ \left( {{w_{a}B} - {w_{b}C}} \right)^{2}} - {\left( {{w_{a}C} - {w_{b}A}} \right)^{2}.}}\end{matrix} & \begin{matrix}(178) \\(179) \\(180) \\\begin{matrix}\; \\(181)\end{matrix}\end{matrix}\end{matrix}$

By defining the constants accordingly, one rewrites as(λ² +R ₁ λ+R ₂)² =R ₃λ²+2R ₄ λ+R ₅.  (182)

Again, one can analyze this equation analogously to the projection ontothe surface maximum slope constraint described above. However,complications arise since there are possibly more than one positive λ's.Instead, one can expand (182)λ⁴+2R ₁λ³+(R ₁ ²+2R ₂ −R ₃)λ²+2(R ₁ R ₂ −R ₄)Δ+R ₂ ² −R ₅=0,  (183)and use the classical Ferrari's method for quartic equation. Then foreach λ>0, defineD=(A−λ)(B+λ)−C ²,  (184)D _(u) =w _(a)(B+λ)−w _(b) C,  (185)D _(v) =w _(b)(A−λ)−w _(a) C.  (186)

Case 2a: D=D_(u)=D_(v)=0. Then one has the systemCu+(B+λ)v=w _(b),  (187)v ² −u ²+1=0.  (188)

Since B>0 and λ>0, the first equation implies v=(w_(b)−Cu)/(B+λ). Thenthe second one becomes:(w _(b) −Cu)²−(B+λ)² u ²+(B+λ)²=0,i.e.,  (189)[C ²−(B+λ)²]u ²−2w _(b) Cu+[w _(b) ²+(B+λ)²]≥0.  (190)If the determinant(w _(b) C)²−(C ²−(B+λ)²)(w _(b) ²+(B+λ)²)≥0,  (191)then two solutions u are obtained. Since there is a unique pair (u, v)with u<0, the quadratic equation (190) must yield one positive and onenegative solutions

$\begin{matrix}{{u = \frac{{w_{b}C} \pm \sqrt{\left( {w_{b}C} \right)^{2} - {\left( {C^{2} - \left( {B + \lambda} \right)^{2}} \right)\left( {w_{b}^{2} + \left( {B + \lambda} \right)^{2}} \right)}}}{C^{2} - \left( {B + \lambda} \right)^{2}}},} & (192)\end{matrix}$and one must also have (C²−(B+λ)²)(w_(b) ²+(B+λ)²)<0. So, the negativesolution u is

$\begin{matrix}{u = {\frac{{w_{b}C} + \sqrt{\left( {w_{b}C} \right)^{2} - {\left( {C^{2} - \left( {B + \lambda} \right)^{2}} \right)\left( {w_{b}^{2} + \left( {B + \lambda} \right)^{2}} \right)}}}{C^{2} - \left( {B + \lambda} \right)^{2}}.}} & (193)\end{matrix}$

Case 2b: D≠0. Then by (175)-(176), u=D_(u)/D and v=D_(v)/D.

Next, among all (u, v)'s, choose the pair with u<0. Then rescalevariables (u, v)←(su, v). Finally, one can use (104)-(106) to obtain(h₁, h₂, h₃).

Multiple Surface Orientation Constraints

It is possible to impose multiple surface orientation constraints on onetriangle. However, it might cause infeasibility. For instance, it iseasy to see that one cannot “tilt” the normal vector toward oppositedirections. Therefore, care must be taken when enforcing multipleorientation constraints. Below is a description of the feasibility insuch cases.

Feasibility Criterion for One Maximum and One Oriented Minimum SlopeConstraints

Given a triangle with unit normal vector {right arrow over (n)}. Supposeconstraints are enforced on {right arrow over (n)} as follows:cos ∠({right arrow over (n)},e ₃)≥1/√{square root over (1+t ²)}(maximumslope)  (194)cos ∠({right arrow over (n)},{right arrow over (q)})≥s/√{square rootover (1+s ²)}(oriented minimum slope)  (195)where e₃=(0,0,1); {right arrow over (q)}=(q₁, q₂, 0) ≡(q₁, q₂) is agiven unit vector; and t, s∈R₊ are also given. Then the problem isfeasible if and only if t≥s.

To verify this, assume {right arrow over (n)}=(u, v, w) is a unit vector(w>0), and its shadow on xy-plane is {right arrow over (n)}₀=(u, v).Then (194) implies

$\begin{matrix}{{\left. {w \geq \frac{1}{\sqrt{1 + t^{2}}}}\Leftrightarrow{u^{2} + v^{2}} \right. = {{{1 - w^{2}} \leq {1 - \frac{1}{1 + t^{2}}}} = \frac{t^{2}}{1 + t^{2}}}},} & (196)\end{matrix}$which is a disk on R² centered at the origin with radius t/√{square rootover (1+t²)}, and (195) implies

$\begin{matrix}{{{{uq}_{1} + {vq}_{2}} \geq \frac{s}{\sqrt{1 + s^{2}}}},} & (197)\end{matrix}$which is a circular segment 2202 with height

$\frac{t}{\sqrt{1 + t^{2}}} - \frac{s}{\sqrt{1 + s^{2}}}$(see FIG. 22). Thus, FIG. 22 illustrates a circular segment defining afeasible set for a maximum slope and an oriented minimum slope inaccordance with one or more embodiments of the invention.

So {right arrow over (n)}₀=(u, v) must lie in the area 2202, which isnonempty if and only if

$\left. {\frac{t}{\sqrt{1 + t^{2}}} \geq \frac{s}{\sqrt{1 + s^{2}}}}\Leftrightarrow{\frac{t^{2}}{1 + t^{2}} \geq \frac{s^{2}}{1 + s^{2}}}\Leftrightarrow{t \geq {s.}} \right.$

Feasibility Criterion for One Maximum and Two Oriented Minimum SlopeConstraints

Given a triangle with unit normal vector {right arrow over (n)}. Suppose{right arrow over (n)} satisfies:

$\begin{matrix}{{{{\cos\;{\angle\left( {\overset{\rightarrow}{n},e_{3}} \right)}} \geq \frac{1}{\sqrt{1 + t^{2}}}},{{\cos\;{\angle\left( {\overset{\rightarrow}{n},{\overset{\rightarrow}{q}}_{1}} \right)}} \geq \frac{s_{1}}{\sqrt{1 + s_{1}^{2}}}},{and}}{{\cos\;{\angle\left( {\overset{\rightarrow}{n},{\overset{\rightarrow}{q}}_{2}} \right)}} \geq {\frac{s_{2}}{\sqrt{1 + s_{2}^{2}}}.}}} & (198)\end{matrix}$where e₃=(0,0,1), {right arrow over (q)}_(i)=(q_(i1), q_(i2), 0)≡(q_(i1)q_(i2)) are unit vectors, and t, s₁, s₂ ∈R₊. Then the problem isfeasible if and only if

$\quad\begin{matrix}{{\left( {{\overset{\rightarrow}{q}}_{1},{\overset{\rightarrow}{q}}_{2}} \right)} \leq {{\cos^{- 1}\left( \frac{s_{1}\sqrt{1 + t^{2}}}{t\sqrt{1 + s_{1}^{2}}} \right)} + {{\cos^{- 1}\left( \frac{s_{2}\sqrt{1 + t^{2}}}{t\sqrt{1 + s_{2}^{2}}} \right)}.}}} & (199)\end{matrix}$

To see this, the (feasible) disk and circular segments (i.e., regions2302 and 2304) are illustrated in FIG. 23, where

${r:=\frac{t}{\sqrt{1 + t^{2}}}},{r_{1}:=\frac{s_{1}}{\sqrt{1 + s_{1}^{2}}}},{{{and}\mspace{14mu} r_{2}}:={\frac{s_{2}}{\sqrt{1 + s_{2}^{2}}}.}}$Thus, the intersection is nonempty if and only if

${\angle\left( {{\overset{->}{q}}_{1},{\overset{->}{q}}_{2}} \right)} \leq {{\cos^{- 1}\left( \frac{r_{1}}{r} \right)} + {{\cos^{- 1}\left( \frac{r_{2}}{r} \right)}.}}$Accordingly, FIG. 23 illustrates a circle representing a feasible regionfor a maximum slope constraint and regions 2302 and 2304 show thefeasible region for two oriented minimum slope constraints.

Feasibility Criterion for Multiple Surface Oriented Minimum SlopeConstraints:

Given a triangle with unit normal vector {right arrow over (n)}, s, t∈R, and a collection of unit vectors Q:={{right arrow over(q)}_(i)}_(i∈1) ⊂R². Let t≥s>0. Suppose the constraints are enforced onn as follows:cos ∠({right arrow over (n)},e ₃)≥1/√{square root over (1+t ²)},  (200)∀i∈I:cos ∠({right arrow over (n)},{right arrow over (q)}_(i))≥s/√{square root over (1+s ²)}.  (201)Then it suffices to enforce 2 constraints in (201) since the rest willbe automatically fulfilled.

To verify this, one can assume {{right arrow over (q)}₁, {right arrowover (q)}₂}⊂Q is the pair that make the largest angle (see FIG. 24) andthe area S is the intersection of two corresponding circular segments.Accordingly, FIG. 24 illustrates multiple surface oriented minimum slopeconstraints via the intersection of two corresponding circular segmentsin accordance with one or more embodiments of the invention.

Similar to previous cases, one can illustrate the disk of radius

$\frac{t}{\sqrt{1 + t^{2}}}$and several circular segments of the same height in FIG. 24

$\left( {{{where}\mspace{14mu} r} = \frac{s}{\sqrt{1 + s^{2}}}} \right).$One can argue that all other {right arrow over (q)}_(i)'s lie in between{right arrow over (q)}₁ and {right arrow over (q)}₂ (otherwise, theintersection with S is empty). Notice that the distances from the originto all circular segments are equal. Thus, one can deduce that Ssatisfies all other oriented minimum slope constraint {right arrow over(q)}_(i). Thus, it suffices to enforce at most two oriented minimumslope constraints.

Further, one can associate each {right arrow over (q)}_(i) withdifferent minimum slope s_(i) and then obtain analogous feasibilitycriterion.

Curvature Minimization

In some design problems, the designer may wish to construct a surfacethat is as “smooth” as possible. This problem is referred to asminimizing the curvature between adjacent triangles in the mesh.

Given the points P_(i)=(p_(i1), p_(i2), h_(i)), i=1, 2, 3, 4, so thatthey form two adjacent triangles Δ₁=P₁P₂P₄ and Δ₂=P₂P₃P₄ (see FIG. 25).In this regard, FIG. 25 illustrates two adjacent triangles with surfacenormal vectors that are used to minimize the curvature in accordancewith one or more embodiments of the invention.

Define a_(i):=p_(i1)−p₄₁ and b_(i):=p_(i2)−p₄₂, for i=1, 2, 3. Then therespective normal vectors are computed by the cross products

$\begin{matrix}{{\overset{\rightarrow}{n}}_{1} = {{\begin{bmatrix}a_{1} \\b_{1} \\{h_{1} - h_{4}}\end{bmatrix} \times \begin{bmatrix}a_{2} \\b_{2} \\{h_{2} - h_{4}}\end{bmatrix}} = \begin{bmatrix}{{{- b_{2}}h_{1}} + {b_{1}h_{2}} + {\left( {b_{2} - b_{1}} \right)h_{4}}} \\{{a_{2}h_{1}} - {a_{1}h_{2}} + {\left( {a_{1} - a_{1}} \right)h_{4}}} \\{{a_{1}b_{2}} - {a_{2}b_{1}}}\end{bmatrix}}} & (202) \\{{\overset{\rightarrow}{n}}_{2} = {{\begin{bmatrix}a_{2} \\b_{2} \\{h_{2} - h_{4}}\end{bmatrix} \times \begin{bmatrix}a_{3} \\b_{3} \\{h_{3} - h_{4}}\end{bmatrix}} = \begin{bmatrix}{{{- b_{3}}h_{2}} + {b_{2}h_{3}} + {\left( {b_{3} - b_{2}} \right)h_{4}}} \\{{a_{3}h_{2}} - {a_{2}h_{3}} + {\left( {a_{2} - a_{3}} \right)h_{4}}} \\{{a_{2}b_{3}} - {a_{3}b_{2}}}\end{bmatrix}}} & (203)\end{matrix}$One can rescale and obtain

$\begin{matrix}{{{\overset{\rightarrow}{n}}_{1} = \left( {\frac{{{- b_{2}}h_{1}} + {b_{1}h_{2}} + {\left( {b_{2} - b_{1}} \right)h_{4}}}{{a_{1}b_{2}} - {a_{2}b_{1}}},\frac{{a_{2}h_{1}} - {a_{1}h_{2}} + {\left( {a_{1} - a_{2}} \right)h_{4}}}{{a_{1}b_{2}} - {a_{2}b_{1}}},1} \right)},} & (204) \\{{\overset{\rightarrow}{n}}_{2} = {\left( {\frac{{{- b_{3}}h_{2}} + {b_{1}h_{3}} + {\left( {b_{3} - b_{2}} \right)h_{4}}}{{a_{2}b_{3}} - {a_{3}b_{2}}},\frac{{a_{3}h_{2}} - {a_{2}h_{3}} + {\left( {a_{2} - a_{3}} \right)h_{4}}}{{a_{2}b_{3}} - {a_{3}b_{2}}},1} \right).}} & (205)\end{matrix}$The curvature can be represented by the difference between {right arrowover (n)}₁ and {right arrow over (n)}₂, i.e.,

$\begin{matrix}{{\overset{\rightarrow}{\delta}}_{12}:={{\overset{\rightarrow}{\delta}}_{\Delta_{1}\Delta_{2}} = {{{\overset{\rightarrow}{n}}_{1} - {\overset{\rightarrow}{n}}_{2}} = {\text{:}\begin{bmatrix}\left\langle {u,h} \right\rangle \\\left\langle {v,h} \right\rangle\end{bmatrix}}}}} & (206)\end{matrix}$where

$\begin{matrix}{{u:=\left( {u_{1},u_{2},u_{3},u_{4}} \right)},{v:=\left( {v_{1},v_{2},v_{3},v_{4}} \right)},{h:=\left( {h_{1},h_{2},h_{3},h_{4}} \right)},} & (207) \\{\mspace{79mu}{{u_{1} = \frac{- b_{2}}{{a_{1}b_{2}} - {a_{2}b_{1}}}},{u_{2} = {\frac{b_{1}}{{a_{1}b_{2}} - {a_{2}b_{1}}} + \frac{b_{3}}{{a_{2}b_{3}} - {a_{3}b_{2}}}}},}} & (208) \\{\mspace{79mu}{{u_{3} = \frac{- b_{2}}{{a_{2}b_{3}} - {a_{3}b_{2}}}},{u_{4} = {- \left( {u_{1} + u_{2} + u_{3}} \right)}},}} & (209) \\{\mspace{79mu}{{v_{1} = \frac{a_{2}}{{a_{1}b_{2}} - {a_{2}b_{1}}}},{v_{2} = {\frac{- a_{1}}{{a_{1}b_{2}} - {a_{2}b_{1}}} + \frac{- a_{3}}{{a_{2}b_{3}} - {a_{3}b_{2}}}}},}} & (210) \\{\mspace{79mu}{{v_{3} = \frac{a_{2}}{{a_{2}b_{3}} - {a_{3}b_{2}}}},{v_{4} = {- {\left( {v_{1} + v_{2} + v_{3}} \right).}}}}} & (211)\end{matrix}$

Now let

be the set of all pairs of adjacent triangles (Δ_(i), Δ_(j)) for whichit is desirable to minimize the curvature. Then for each pair (Δ_(i),Δ_(j))∈

, one finds the corresponding vectors h_(ij)=(h₁ ^(ij), h₂ ^(ij), h₃^(ij), h₄ ^(ij))∈R⁴, u_(ij), v_(ij)∈R⁴, and computes the correspondingcurvature δ_(ij)=(

u_(ij), h_(ij)

,

v_(ij), h_(ij)

). One can then aim to minimize all the computed curvatures. Thus, onecan arrive at the objective

$\begin{matrix}{{G_{1,*}(x)}:={\sum\limits_{{({\Delta_{i},\Delta_{j}})} \in \mathcal{Z}}\;{{\overset{\rightarrow}{\delta}}_{ij}}_{*}}} & (212)\end{matrix}$where ∥⋅∥_(*) can be either 1-norm or max-norm in R². For 1-norm, theobjective is

$\begin{matrix}{{{G_{1,1}(x)}:={{\sum\limits_{{({\Delta_{i},\Delta_{j}})} \in \mathcal{Z}}\;{\left\langle {u_{ij},h_{ij}} \right\rangle }} + {\left\langle {v_{ij},h_{ij}} \right\rangle }}},} & (213)\end{matrix}$and for max-norm, the objective is

$\begin{matrix}{{G_{1,\infty}(x)}:={\sum\limits_{{({\Delta_{i},\Delta_{j}})} \in \mathcal{Z}}{\max{\left\{ {{\left\langle {u_{ij},h_{ij}} \right\rangle } + {\left\langle {v_{ij},h_{ij}} \right\rangle }} \right\}.}}}} & (214)\end{matrix}$

Simplified Computations for Symmetric Cases.

Suppose the two dimensional mesh satisfies the following symmetry: forevery pair of adjacent triangles (P₁P₂P₄, P₂P₃P₄)∈

, there exists t∈R such that{right arrow over (P ₄ P ₁)}+{right arrow over (P ₄ P ₃)}=t{right arrowover (P ₄ P ₂)}.  (215)Then it follows that (a₁, b₁₎+(a₃, b₃)=t(a₂, b₂). So one can deducea₁b₂−a₂b₁=a₂b₃−a₃b₂. Using (207)-(211), one obtains

$\begin{matrix}{u = {{\frac{- b_{2}}{{a_{1}b_{2}} - {a_{2}b_{1}}}\left( {1,{- t},1,{t - 2}} \right)\mspace{14mu}{and}\mspace{14mu} v} = {\frac{a_{2}}{{a_{1}b_{2}} - {a_{2}b_{1}}}{\left( {1,{- t},1,{t - 2}} \right).}}}} & (216)\end{matrix}$That means u and v are parallel. This simplifies the computations for(213) and (214).

Proximity Operators

Since optimization methods can require proximity operators, one willderive the necessary formulas. It is sometimes convenient to compute theproximity operator P_(ƒ) via the proximity operator of its Fenchelconjugate ƒ*, which is defined by ƒ*:X→R:x*

sup_(x∈X)(

x*, x

−ƒ(x)). Indeed, if γ>0, then (see, e.g., [2,

Theorem 14.3(ii)])∀x∈X:x=P _(γƒ)(x)+γP _(γ) ⁻¹ _(ƒ*(γ) ⁻¹ x).  (217)One may also recall a useful formula from [2, Lemma 2.3]: if ƒ:X→R isconvex and positively homogeneous, α>0, w∈X, and h:X→R:x

αƒ(x−w), then

$\begin{matrix}{{P_{h}(x)} = {{w + {\alpha\;{P_{f}\left( \frac{x - w}{\alpha} \right)}}} = {x - {\alpha\;{{P_{f^{*}}\left( \frac{x - w}{\alpha} \right)}.}}}}} & (218)\end{matrix}$

Formula for Proximity Operators:

Let {u_(i)}_(i∈I) be a system of finitely many vectors in R^(n), and

$\begin{matrix}\left. {f\text{:}\mspace{14mu} R^{''}}\rightarrow\left. {R\text{:}\mspace{14mu} x}\rightarrow{\max\limits_{i \in I}{\left\{ {\left\langle {u_{i},x} \right\rangle } \right\}.}} \right. \right. & (219)\end{matrix}$

Then ƒ*=t_(D) where D:=convU_(i∈I){u_(i),−u_(i)}. Consequently, theproximity operator of ƒ can be computed by P_(ƒ)=Id−P_(D), where P_(D)is the projection onto the set D.

To see this, first, suppose x*∈D, then one can express

$\begin{matrix}{x^{*} = {{\sum\limits_{i \in I}\;{\lambda_{i}u_{i}\mspace{14mu}{where}\mspace{14mu}{\sum\limits_{i \in I}\;{\lambda_{i}}}}} \leq 1.}} & (220)\end{matrix}$It follows that for all x∈X,

$\begin{matrix}{{\left\langle {x^{*},x} \right\rangle - {f(x)}} = {{{\sum\limits_{i \in I}{\lambda_{i}\left\langle {u_{i},x} \right\rangle}} - {\max\limits_{i \in I}{\left\langle {u_{i},x} \right\rangle }}} \leq {\left( {{\sum\limits_{i \in I}\lambda_{i}} - 1} \right){\max\limits_{i \in I}{\left\langle {u_{i},x} \right\rangle }}} \leq 0.}} & (221)\end{matrix}$So, ƒ* (x*)=sup_(x∈X)[

x*, x

−ƒ(x)]≤0. Notice that equality happens if one sets x=0. Thus, ƒ*(x*)=0.

Now suppose x*∉D. Since D is nonempty, closed, and convex, the classicseparation theorem implies that there exists x∈X such that

x*,x

>

u,x

for all u∈D  (222)This leads to

x*, x

−ƒ(x)>0. Since ƒ is homogeneous, one can haveƒ*(x*)≥

x*,λx

−ƒ(λx)→+∞ as λ→+∞.  (223)So ƒ*(x*)=+∞. Therefore, one concludes that ƒ*=t_(D). Now employingformula (218) and noting that P_(ƒ)*=P_(D), the projection onto D (seethe beginning of the iterative projection methods optimization problemsdescription above), one obtains the proximity operator for ƒ.

Next, one presents two special cases of the above proximity formula,which will be used in curvature minimization.

Example 1

Given α>0, a vector u ∈R^(n) \{0} and the functionƒ:R ^(n) →R: x

α|

u,x

|.  (224)Then the proximity operator of ƒ isP _(ƒ) =Id−P _(D) where D:=[−αu,αu].  (225)The explicit projection onto D is given below (see also, [4, Theorem2.7]).

$\begin{matrix}{{P_{D}x} = {{\min\left\{ {1,{\max\left\{ {{- 1},\frac{\left\langle {{\alpha\; u},x} \right\rangle}{{{\alpha\; u}}^{2}}} \right\}}} \right\}\alpha\; u} = {\min\left\{ {\alpha,{\max\left\{ {{- \alpha},\frac{\left\langle {u,x} \right\rangle}{{u}^{2}}} \right\}}} \right\}{u.}}}} & (226)\end{matrix}$

Example 2

Given two vectors u, v∈R^(n) andƒ:R ^(n) →R:x

max{|

u,x

|,|

v,x

|}.  (227)Then the proximity operator of ƒ isP _(ƒ) =Id−P _(D) where D:=conv{u,v,−u,−v}.  (228)In general, P_(D) is the projection onto a parallelogram in R^(n).Approximation of Surface Slope Constraint Projections

In this section, the problem of surface slope constraint projections areconsidered. Such projections may be written in the form

$\begin{matrix}{{\min\limits_{u \in R^{2}}{\varphi(u)}}:={{\frac{A}{2}u_{1}^{2}} + {{Cu}_{1}u_{2}} + {\frac{B}{2}u_{2}^{2}} - {w_{a}u_{1}} - {w_{b}u_{2}}}} & (229)\end{matrix}$subject to u∈Ω:={u∈R ² |g(u)≤0}.  (230)

If one defines

${M:={{\begin{bmatrix}A & C \\C & B\end{bmatrix}\mspace{14mu}{and}\mspace{14mu} w}:=\begin{bmatrix}W_{a} \\W_{b}\end{bmatrix}}},$then

$\begin{matrix}{{\varphi(u)} = {{\frac{1}{2}u^{T}{Mu}} - {w^{T}{u.}}}} & (231)\end{matrix}$In practice, solving (229)-(230) directly as described above (e.g., asdescribed above in the “Projection onto the surface maximum slopeconstraint” section and the “Projection onto Surface Oriented MinimumSlope Constraint” section) might not be straightforward as one mayencounter numerical instability. Therefore, one can use an alternativeapproach where the constraint Ω is approximated by a regular polygon.

Approximation for Maximum Slope Constraint

Consider the projection onto maximum slope constraint written in theform

$\begin{matrix}{{\min\limits_{u \in R^{2}}{\varphi(u)}}:={{\frac{1}{2}u^{T}{Mu}} - {w^{T}u}}} & (232) \\{{{{subject}\mspace{14mu}{to}\mspace{14mu} u} \in \Omega}:=\left\{ {\left( {u_{1},u_{2}} \right) \in {R^{2}{\left. {{u_{1}^{2} + u_{2}^{2}} \leq 1} \right\}.}}} \right.} & (233)\end{matrix}$The unit disk Ω may be replaced by a convex (or simple) regular polygonthat is contained in Ω. For example, the unit disk may be replaced by ahexagon (see FIG. 26). In this regard, FIG. 26 illustrates a hexagonthat is used to replace a unit disk for a maximum slope constraint inaccordance with one or more embodiments of the invention. Morespecifically, let m≥3 be an integer and 0≤θ₁≤θ₂≤ . . . ≤θ_(m)<2π. Definethe points on the unit circle p_(k):=(cos θ_(k), sin θ_(k)) for k=1, 2,. . . , m, and set also p₀:=p_(m). Now define the convex m-polygonU:=conv{p_(k)|k=0, 1, . . . , m}. One then arrives at the approximateproblem

$\begin{matrix}{{\min\limits_{u \in R^{2}}{\varphi(u)}}:={{{\frac{1}{2}u^{T}{Mu}} - {w^{T}u\mspace{14mu}{subject}\mspace{14mu}{to}\mspace{14mu} u}} \in {U.}}} & (234)\end{matrix}$This problem may be solved as follows. Let u₀:=M⁻¹w be the vertex of theparaboloid φ(u).

Case 1: u₀ ∈Ω. Then u₀ is even the solution of (232)-(233), thus, theproblem is solved.

Case 2: u₀ ∉C. So one also has u₀∉U. The solution of (234) may be foundby the following steps.

Step 2a: Find the index setK:={k|0 and u ₀ are on opposite sides of the line p _(k) p_(k+1)}.  (235)

Step 2b: Minimize φ(u) on each edge p_(k)p_(k+1) for k∈K. Consider φ(u)on an edge p_(k)p_(k+1). Define Δp_(k)=p_(k+1)−p_(k). Take a pointu∈p_(k)p_(k+1), thenu=p _(k) +tΔp _(k) for some t∈[0,1].  (236)So

$\begin{matrix}{{\psi(t)}:={{\varphi(u)} = {{\frac{1}{2}\left( {p_{k} + {t\;\Delta\; p_{k}}} \right)^{T}{M\left( {p_{k} + {t\;\Delta\; p_{k}}} \right)}} - {w^{T}\left( {p_{k} + {t\;\Delta\; p_{k}}} \right)}}}} & {{~~~~~}(237)} \\{= {{\frac{t^{2}}{2}\Delta\; p_{k}^{T}M\;\Delta\; p_{k}} - {t\left( {{w^{T}\Delta\; p_{k}} - {p_{k}^{T}M\;\Delta\; p_{k}}} \right)} + {\frac{1}{2}p_{k}^{T}{Mp}_{k}} -}} & {(238)} \\{w^{T}{p_{k}.}} & \end{matrix}$The vertex of ψ(t) is

${\overset{\_}{t}:=\frac{{w^{T}\Delta\; p_{k}} - {p_{k}^{T}M\;\Delta\; p_{k}}}{\Delta\; p_{k}^{T}M\;\Delta\; p_{k}}},$so the minimizer of ψ(t) on [0,1] ist _(k):=min{1,max{0,t}}.  (239)

Step 2c: If there exists k∈K such that 0<t_(k)<1, thenu=p_(k)+t_(k)Δp_(k) is the solution of (234). Otherwise, the solution of(234) is

$\begin{matrix}{u = {{p_{k_{0}} + {t_{k_{0}}\Delta\; p_{k_{0}}\mspace{14mu}{where}\mspace{14mu} k_{0}}}:={\underset{k \in K}{argmin}{{\varphi\left( {p_{k} + {t_{k}\Delta\; p_{k}}} \right)}.}}}} & (240)\end{matrix}$

Hence, the projection to the surface maximum slope constraint isapproximated by using a regular polygon instead of a unit disk.

Approximation for Oriented Minimum Slope Constraint

Consider the projection onto oriented minimum slope constraint writtenin the form

$\begin{matrix}{\mspace{79mu}{{\min\limits_{u \in R^{2}}{\varphi(u)}}:={{\frac{1}{2}u^{T}{Mu}} - {w^{T}u}}}} & (241) \\{{{{subject}\mspace{14mu}{to}\mspace{14mu} u} \in \Omega}:={\left\{ {\left. {\left( {u_{1},u_{2}} \right) \in R^{2}} \middle| {u_{1} < 0} \right.,{{u_{2}^{2} - \frac{u_{1}^{2}}{s^{2}} + 1} \leq 0}} \right\}.}} & (242)\end{matrix}$The constraint Ω may be replaced byU:={(u ₁ ,u ₂)∈R ² |u ₁ ≤−s(|u ₂|+1)}.  (243)FIG. 27 shows how the original half cone is now replaced with theabsolute value function in accordance with one or more embodiments ofthe invention. One can easily check that U⊂Ω. So one arrives at anapproximate problem

$\begin{matrix}{{\min\limits_{u \in R^{2}}{\varphi(u)}}:={{{\frac{1}{2}u^{T}{Mu}} - {w^{T}u\mspace{14mu}{subject}\mspace{14mu}{to}\mspace{14mu} u}} \in {U.}}} & (244)\end{matrix}$Let u₀:=M⁻¹w be the vertex of the paraboloid φ(u).

Case 1: u₀ ∈C. Then u₀ is even the solution of (241)-(242), thus, theproblem is solved.

Case 2: u₀ ∈C. So u₀ ∈U. One solves (244) by first finding the minimizerof φ(u) on the boundary of U, which consists of two half-lines{p+tΔp ₁ |t∈R ₊} and {p+tΔp ₂ |tΣR ₊},  (245)where p:=(−s, 0), Δp₁:=(−s, 1), and Δp₂:=(−s, −1).Let u:=p+tΔp for t∈[0, +∞). Similar to the computation above, one has

$\begin{matrix}{{\psi(t)}:={{\varphi(u)} = {{\frac{t^{2}}{2}\Delta\; p^{T}M\;\Delta\; p} - {t\left( {{w^{T}\Delta\; p} - {p^{T}M\;\Delta\; p}} \right)} + {\frac{1}{2}p^{T}{Mp}} - {w^{T}{p.}}}}} & (246)\end{matrix}$So the minimizer of φ(u) on [0, +∞) is

$\begin{matrix}{t:={\max{\left\{ {0,\frac{{w^{T}\Delta\; p} - {p^{T}M\;\Delta\; p}}{\Delta\; p^{T}M\;\Delta\; p}} \right\}.}}} & (247)\end{matrix}$Second, one sets Δp=Δp₁ and Δp=Δp₂ in the above formula to obtain t₁ andt₂, respectively. Finally, comparing the values φ(p+t₁Δp₁) andφ(p+t₂Δp₂), one obtains the minimizer of (244).

Visualization Examples

The closed-form proximity and projection operators can be used in the DRmethod to solve various design problems described above.

FIG. 28 shows the starting situation for the problem shown in FIG. 7 inaccordance with one or more embodiments of the invention. Triangles 2802represent triangles that violate the design constraints, and triangles2902 (from FIG. 29) represent triangles that satisfy the designconstraints. FIG. 29 shows the mesh after 10 iterations in accordancewith one or more embodiments of the invention. FIG. 30 shows the meshafter 29 iterations, and FIG. 31 shows the final mesh after convergencein accordance with one or more embodiments of the invention.

Consider FIG. 28 again, this time as a starting situation for theproblem shown in FIG. 8. FIG. 32 is after 15 iterations, FIG. 33 isafter 100 iterations, and FIG. 34 is the final solution.

FIG. 35 is the starting situation for the problem shown in FIG. 9. FIG.36 is after 10 iterations, FIG. 37 is after 20 iterations, and FIG. 38is the final solution.

Logical Flow

FIG. 39 illustrates the logical flow for designing a surface inaccordance with one or more embodiments of the invention.

At step 3902, a triangular surface mesh representative of an existingsurface is obtained. The triangular surface mesh includes two or moretriangles that are connected by vertices and edges. In embodiments ofthe invention, the surface may be a land surface (e.g., that is actuallyconstructed based on the design [e.g., a computer-aided design—CAD]).

At step 3904, one or more design constraint sets are determined based onone or more design constraints. The design constraints can includevarious constraints including a maximum slope constraint for a firsttriangle of the two or more triangles in the triangular surface mesh.The maximum slope constraint is a maximum angle between a normal vectorof the first triangle and a reference vector (e.g., a Z-vectorrepresentative of gravity).

An additional design constraint may be a surface oriented minimum slopeconstraint for the first triangle. The surface oriented minimum slopeconstraint is a minimum angle between the normal vector of the firsttriangle and a directional vector.

At step 3906, one or more heights of the vertices of the first triangleare projected onto the one or more design constraint sets such that thenormal vector satisfies all of the one or more design constraints. Theprojecting includes modifying the one or more heights by a minimumEuclidian distance. Further, the projecting may be performed iterativelyuntil all of the one or more design constraints are satisfied. Inaddition, the projecting onto the design constraint set for the maximumslope constraint may include moving the normal vector onto a conedefined by an origin of the reference vector on the triangular surfacemesh and a unit circle around the reference vector that is defined bythe maximum angle. Alternatively, or in addition, the projecting onto tothe design constraint set for the maximum slope constraint may includemoving the normal vector onto a polygonal shaped pyramid defined by thereference vector on the triangular surface mesh and a regular polygonaround the reference vector, wherein edges of the polygonal shapedpyramid are at the maximum angle with the reference vector.

At step 3908, a design of the surface represented by the triangularsurface mesh is generated (based on the projecting).

Hardware Environment

FIG. 40 is an exemplary hardware and software environment 4000 used toimplement one or more embodiments of the invention. The hardware andsoftware environment includes a computer 4002 and may includeperipherals. Computer 4002 may be a user/client computer, servercomputer, or may be a database computer. The computer 4002 comprises ageneral purpose hardware processor 4004A and/or a special purposehardware processor 4004B (hereinafter alternatively collectivelyreferred to as processor 4004) and a memory 4006, such as random accessmemory (RAM). The computer 4002 may be coupled to, and/or integratedwith, other devices, including input/output (I/O) devices such as akeyboard 4014, a cursor control device 4016 (e.g., a mouse, a pointingdevice, pen and tablet, touch screen, multi-touch device, etc.) and aprinter 4028. In one or more embodiments, computer 4002 may be coupledto, or may comprise, a portable or media viewing/listening device 4032(e.g., an MP3 player, IPOD, NOOK, portable digital video player,cellular device, personal digital assistant, etc.). In yet anotherembodiment, the computer 4002 may comprise a multi-touch device, mobilephone, gaming system, internet enabled television, television set topbox, or other internet enabled device executing on various platforms andoperating systems.

In one embodiment, the computer 4002 operates by the general purposeprocessor 4004A performing instructions defined by the computer program4010 under control of an operating system 4008. The computer program4010 and/or the operating system 4008 may be stored in the memory 4006and may interface with the user and/or other devices to accept input andcommands and, based on such input and commands and the instructionsdefined by the computer program 4010 and operating system 4008, toprovide output and results.

Output/results may be presented on the display 4022 or provided toanother device for presentation or further processing or action. In oneembodiment, the display 4022 comprises a liquid crystal display (LCD)having a plurality of separately addressable liquid crystals.Alternatively, the display 4022 may comprise a light emitting diode(LED) display having clusters of red, green and blue diodes driventogether to form full-color pixels. Each liquid crystal or pixel of thedisplay 4022 changes to an opaque or translucent state to form a part ofthe image on the display in response to the data or informationgenerated by the processor 4004 from the application of the instructionsof the computer program 4010 and/or operating system 4008 to the inputand commands. The image may be provided through a graphical userinterface (GUI) module 4018. Although the GUI module 4018 is depicted asa separate module, the instructions performing the GUI functions can beresident or distributed in the operating system 4008, the computerprogram 4010, or implemented with special purpose memory and processors.

In one or more embodiments, the display 4022 is integrated with/into thecomputer 4002 and comprises a multi-touch device having a touch sensingsurface (e.g., track pod or touch screen) with the ability to recognizethe presence of two or more points of contact with the surface. Examplesof multi-touch devices include mobile devices (e.g., IPHONE, NEXUS S,DROID devices, etc.), tablet computers (e.g., IPAD, HP TOUCHPAD),portable/handheld game/music/video player/console devices (e.g., IPODTOUCH, MP3 players, NINTENDO 3DS, PLAYSTATION PORTABLE, etc.), touchtables, and walls (e.g., where an image is projected through acrylicand/or glass, and the image is then backlit with LEDs).

Some or all of the operations performed by the computer 4002 accordingto the computer program 4010 instructions may be implemented in aspecial purpose processor 4004B. In this embodiment, some or all of thecomputer program 4010 instructions may be implemented via firmwareinstructions stored in a read only memory (ROM), a programmable readonly memory (PROM) or flash memory within the special purpose processor4004B or in memory 4006. The special purpose processor 4004B may also behardwired through circuit design to perform some or all of theoperations to implement the present invention. Further, the specialpurpose processor 4004B may be a hybrid processor, which includesdedicated circuitry for performing a subset of functions, and othercircuits for performing more general functions such as responding tocomputer program 4010 instructions. In one embodiment, the specialpurpose processor 4004B is an application specific integrated circuit(ASIC).

The computer 4002 may also implement a compiler 4012 that allows anapplication or computer program 4010 written in a programming languagesuch as C, C++, Assembly, SQL, PYTHON, PROLOG, MATLAB, RUBY, RAILS,HASKELL, or other language to be translated into processor 4004 readablecode. Alternatively, the compiler 4012 may be an interpreter thatexecutes instructions/source code directly, translates source code intoan intermediate representation that is executed, or that executes storedprecompiled code. Such source code may be written in a variety ofprogramming languages such as JAVA, JAVASCRIPT, PERL, BASIC, etc. Aftercompletion, the application or computer program 4010 accesses andmanipulates data accepted from I/O devices and stored in the memory 4006of the computer 4002 using the relationships and logic that weregenerated using the compiler 4012.

The computer 4002 also optionally comprises an external communicationdevice such as a modem, satellite link, Ethernet card, or other devicefor accepting input from, and providing output to, other computers 4002.

In one embodiment, instructions implementing the operating system 4008,the computer program 4010, and the compiler 4012 are tangibly embodiedin a non-transitory computer-readable medium, e.g., data storage device4020, which could include one or more fixed or removable data storagedevices, such as a zip drive, floppy disc drive 4024, hard drive, CD-ROMdrive, tape drive, etc. Further, the operating system 4008 and thecomputer program 4010 are comprised of computer program 4010instructions which, when accessed, read and executed by the computer4002, cause the computer 4002 to perform the steps necessary toimplement and/or use the present invention or to load the program ofinstructions into a memory 4006, thus creating a special purpose datastructure causing the computer 4002 to operate as a specially programmedcomputer executing the method steps described herein. Computer program4010 and/or operating instructions may also be tangibly embodied inmemory 4006 and/or data communications devices 4030, thereby making acomputer program product or article of manufacture according to theinvention. As such, the terms “article of manufacture,” “program storagedevice,” and “computer program product,” as used herein, are intended toencompass a computer program accessible from any computer readabledevice or media.

Of course, those skilled in the art will recognize that any combinationof the above components, or any number of different components,peripherals, and other devices, may be used with the computer 4002.

FIG. 41 schematically illustrates a typical distributed/cloud-basedcomputer system 4100 using a network 4104 to connect client computers4102 to server computers 4106. A typical combination of resources mayinclude a network 4104 comprising the Internet, LANs (local areanetworks), WANs (wide area networks), SNA (systems network architecture)networks, or the like, clients 4102 that are personal computers orworkstations (as set forth in FIG. 40), and servers 4106 that arepersonal computers, workstations, minicomputers, or mainframes (as setforth in FIG. 40). However, it may be noted that different networks suchas a cellular network (e.g., GSM [global system for mobilecommunications] or otherwise), a satellite based network, or any othertype of network may be used to connect clients 4102 and servers 4106 inaccordance with embodiments of the invention.

A network 4104 such as the Internet connects clients 4102 to servercomputers 4106. Network 4104 may utilize ethernet, coaxial cable,wireless communications, radio frequency (RF), etc. to connect andprovide the communication between clients 4102 and servers 4106.Further, in a cloud-based computing system, resources (e.g., storage,processors, applications, memory, infrastructure, etc.) in clients 4102and server computers 4106 may be shared by clients 4102, servercomputers 4106, and users across one or more networks. Resources may beshared by multiple users and can be dynamically reallocated per demand.In this regard, cloud computing may be referred to as a model forenabling access to a shared pool of configurable computing resources.

Clients 4102 may execute a client application or web browser andcommunicate with server computers 4106 executing web servers 4110. Sucha web browser is typically a program such as MICROSOFT INTERNETEXPLORER, MOZILLA FIREFOX, OPERA, APPLE SAFARI, GOOGLE CHROME, etc.Further, the software executing on clients 4102 may be downloaded fromserver computer 4106 to client computers 4102 and installed as a plug-inor ACTIVEX control of a web browser. Accordingly, clients 4102 mayutilize ACTIVEX components/component object model (COM) or distributedCOM (DCOM) components to provide a user interface on a display of client4102. The web server 4110 is typically a program such as MICROSOFT'SINTERNET INFORMATION SERVER.

Web server 4110 may host an Active Server Page (ASP) or Internet ServerApplication Programming Interface (ISAPI) application 4112, which may beexecuting scripts. The scripts invoke objects that execute businesslogic (referred to as business objects). The business objects thenmanipulate data in database 4116 through a database management system(DBMS) 4114. Alternatively, database 4116 may be part of, or connecteddirectly to, client 4102 instead of communicating/obtaining theinformation from database 4116 across network 4104. When a developerencapsulates the business functionality into objects, the system may bereferred to as a component object model (COM) system. Accordingly, thescripts executing on web server 4110 (and/or application 4112) invokeCOM objects that implement the business logic. Further, server 4106 mayutilize MICROSOFT'S TRANSACTION SERVER (MTS) to access required datastored in database 4116 via an interface such as ADO (Active DataObjects), OLE DB (Object Linking and Embedding DataBase), or ODBC (OpenDataBase Connectivity).

Generally, these components 4100-4116 all comprise logic and/or datathat is embodied in/or retrievable from device, medium, signal, orcarrier, e.g., a data storage device, a data communications device, aremote computer or device coupled to the computer via a network or viaanother data communications device, etc. Moreover, this logic and/ordata, when read, executed, and/or interpreted, results in the stepsnecessary to implement and/or use the present invention being performed.

Although the terms “user computer”, “client computer”, and/or “servercomputer” are referred to herein, it is understood that such computers4102 and 4106 may be interchangeable and may further include thin clientdevices with limited or full processing capabilities, portable devicessuch as cell phones, notebook computers, pocket computers, multi-touchdevices, and/or any other devices with suitable processing,communication, and input/output capability.

Of course, those skilled in the art will recognize that any combinationof the above components, or any number of different components,peripherals, and other devices, may be used with computers 4102 and4106. Further, embodiments of the invention are implemented as asoftware application on a client 4102 or server computer 4106. Inaddition, as described above, the client 4102 or server computer 4106may comprise a thin client device or a portable device that has amulti-touch-based display.

CONCLUSION

This concludes the description of the preferred embodiment of theinvention. The following describes some alternative embodiments foraccomplishing the present invention. For example, any type of computer,such as a mainframe, minicomputer, or personal computer, or computerconfiguration, such as a timesharing mainframe, local area network, orstandalone personal computer, could be used with the present invention.

The foregoing description of the preferred embodiment of the inventionhas been presented for the purposes of illustration and description. Itis not intended to be exhaustive or to limit the invention to theprecise form disclosed. Many modifications and variations are possiblein light of the above teaching. It is intended that the scope of theinvention be limited not by this detailed description, but rather by theclaims appended hereto.

REFERENCES

-   [1] BAUSCHKE, H. H., AND BORWEIN, J. M. On projection algorithms for    solving convex feasibility problems. SIAM Rev. 38, 3 (1996),    367-426.-   [2] BAUSCHKE, H. H., AND COMBETTES, P. L. Convex analysis and    monotone operator theory in Hilbert spaces. CMS Books in    Mathematics/Ouvrages de Mathématiques de la SMC. Springer, New    York, 2011. With a foreword by Hédy Attouch.-   [3] BAUSCHKE, H. H., AND KOCH, V. R. Projection methods: Swiss army    knives for solving feasibility and best approximation problems with    halfspaces. In Infinite products of operators and their    applications, vol. 636 of Contemp. Math. Amer. Math. Soc.,    Providence, R.I., 2015, pp. 1-40.-   [4] BAUSCHKE, H. H., KOCH, V. R., AND PHAN, H. M. Stadium norm and    Douglas-Rachford splitting: a new approach to road design    optimization. Oper. Res. 64, 1 (2016), 201-218.-   [5] BINIAZ, A., AND DASTGHAIBYFARD, G. Drainage reality in terrains    with higher-order Delaunay triangulations. Springer Berlin    Heidelberg, Berlin, Heidelberg, 2008, pp. 199-211.-   [6] BINIAZ, A., AND DASTGHAIBYFARD, G. Slope fidelity in terrains    with higher-order delaunay triangulations. In 16th International    Conference in Central Europe on Computer Graphics, Visualization,    and Computer Vision (Plzen, Czech Republic, 2008), Václav    Skala-UNION Agency, pp. 17-23.-   [7] BOT, R. I., CSETNEK, E. R., AND HEINRICH, A. A primal-dual    splitting algorithm for finding zeros of sums of maximal monotone    operators. SIAM J. Optim. 23, 4 (2013), 2011-2036.-   [8] BRICEÑO ARIAS, L. M., AND COMBETTES, P. L. A monotone+skew    splitting model for composite monotone inclusions in duality.    SIAM J. Optim. 21, 4 (2011), 1230-1250.-   [9] CARDANO, G. The great art or The rules of algebra. The M.I.T.    Press, Cambridge, Mass.-London, 1968. Translated from the Latin and    edited by T. Richard Witmer, With a foreword by Oystein Ore.-   [10] CEGIELSKI, A. Iterative methods for fixed point problems in    Hilbert spaces, vol. 2057 of Lecture Notes in Mathematics. Springer,    Heidelberg, 2012.-   [11] CENSOR, Y., CHEN, W., COMBETTES, P. L., DAVIDI, R., AND    HERMAN, G. T. On the effectiveness of projection methods for convex    feasibility problems with linear inequality constraints. Comput.    Optim. Appl. 51, 3 (2012), 1065-1088.-   [12] CENSOR, Y., AND ZENIOS, S. A. Parallel optimization. Numerical    Mathematics and Scientific Computation. Oxford University Press, New    York, 1997. Theory, algorithms, and applications, With a foreword by    George B. Dantzig.-   [13] CORNISH, G., AND MUIR GRAVES, R. Classic Golf Hole Design. John    Wiley & Sons, Inc., Hoboken DC, 2002.-   [14] DOUGLAS, JR., J., AND RACHFORD, JR., H. H. On the numerical    solution of heat conduction problems in two and three space    variables. Trans. Amer. Math. Soc. 82 (1956), 421-439.-   [15] HAWTREE, F. W. The Golf Course, Planning design, construction    and maintenance. E & FN Spon, and imprint of Chapman & Hall, London,    U K, 2005.-   [16] IRVING, R. Beyond the Quadratic Formula. The Mathematical    Association of America, Washington D.C., 2010.-   [17] ITOH, T., AND SHIMADA, K. Automatic conversion of triangular    meshes into quadrilateral meshes with directionality. International    Journal of CADCAM 1, 1 (2002), 11-21.-   [18] JOHNSTON, B. P., SULLIVAN, J. M., AND KWASNIK, A. Automatic    conversion of triangular finite element meshes to quadrilateral    elements. International Journal for Numerical Methods in Engineering    31, 1 (1991), 67-84.-   [19] KARUSH, W. MINIMA OF FUNCTIONS OF SEVERAL VARIABLES WITH    INEQUALITIES AS SIDE CONDITIONS. ProQuest LLC, Ann Arbor,    Mich., 1939. Thesis (SM)—The University of Chicago.-   [20] KUHN, H. W., AND TUCKER, A. W. Nonlinear programming. In    Proceedings of the Second Berkeley Symposium on Mathematical    Statistics and Probability, 1950 (1951), University of California    Press, Berkeley and Los Angeles, pp. 481-492.-   [21] LIONS, P.-L., AND MERCIER, B. Splitting algorithms for the sum    of two nonlinear operators. SIAM J. Numer. Anal. 16, 6 (1979),    964-979.-   [22] MOORE, I.D., GRAYSON, R., AND LADSON, A. Digital terrain    modelling: a review of hydrological, geomorphological, and    biological applications. Hydrological processes 5, 1 (1991), 3-30.-   [23] MOREAU, J.-J. Proximité et dualité dans un espace hilbertien.    Bull. Soc. Math. France 93 (1965), 273-299.-   [24] RANK, E., SCHWEINGRUBER, M., AND SOMMER, M. Adaptive mesh    generation and transformation of triangular to quadrilateral meshes.    Communications in Numerical Methods in Engineering 9, 2 (1993),    121-129.-   [25] ROCKAFELLAR, R. T. Convex analysis. Princeton Mathematical    Series, No. 28. Princeton University Press, Princeton, N.J., 1970.-   [26] RUSZCZYŃSKI, A. Nonlinear optimization. Princeton University    Press, Princeton, N.J., 2006.

What is claimed is:
 1. A computer-implemented method for designing asurface, comprising: (a) obtaining, in a computer, a triangular surfacemesh representative of an existing physical surface, wherein thetriangular surface mesh comprises two or more triangles that areconnected by vertices and edges; (b) determining, in the computer, oneor more design constraint sets based on one or more design constraints,wherein at least one of the one or more design constraints comprise: amaximum slope constraint for a first triangle of the two or moretriangles in the triangular surface mesh, wherein the maximum slopeconstraint comprises a maximum angle between a normal vector of thefirst triangle and a reference vector; (c) projecting, in the computer,one or more heights of the vertices of the first triangle onto the oneor more design constraint sets such that the normal vector satisfies allof the one or more design constraints, wherein the projecting comprisesmodifying the one or more heights by a minimum Euclidian distance, andwherein the modifying changes the triangular surface mesh into aconverged triangular surface mesh; (d) generating, in the computer, adesign of the surface represented by the converged triangular surfacemesh, wherein the generating is based on the projecting; and (e)rendering the design of the surface in a computer-aided design (CAD)application, wherein the surface is constructed based on the design. 2.The computer-implemented method of claim 1, wherein the one or moredesign constraints further comprise: a surface oriented minimum slopeconstraint for the first triangle, wherein the surface oriented minimumslope constraint comprises a minimum angle between the normal vector ofthe first triangle and a directional vector.
 3. The computer-implementedmethod of claim 1, wherein: the surface comprises a land surface; andthe method further comprises constructing the land surface based on thedesign.
 4. The computer-implemented method of claim 3, wherein: thedesign comprises a computer-aided design (CAD) for the land surface. 5.The computer-implemented method of claim 1, wherein the reference vectorcomprises a Z-vector representative of gravity.
 6. Thecomputer-implemented method of claim 1, wherein: the projecting isperformed iteratively until all of the one or more design constraintsare satisfied.
 7. The computer-implemented method of claim 1, wherein:the projecting onto the design constraint set for the maximum slopeconstraint comprises moving the normal vector onto a cone defined by anorigin of the reference vector on the triangular surface mesh and a unitcircle around the reference vector that is defined by the maximum angle.8. The computer-implemented method of claim 1, wherein: the projectingonto to the design constraint set for the maximum slope constraintcomprises moving the normal vector onto a polygonal shaped pyramiddefined by the reference vector on the triangular surface mesh and aregular polygon around the reference vector, wherein edges of thepolygonal shaped pyramid are at the maximum angle with the referencevector.
 9. A system for designing a surface, comprising: (a) a computerhaving a memory and a processor; (b) a computer aided design (CAD)application executed by the processor; (c) a triangular surface mesh,obtained in the CAD application, wherein the triangular surface mesh isrepresentative of an existing physical surface and comprises two or moretriangles that are connected by vertices and edges; (d) one or moredesign constraint sets, determined by the CAD application, based on oneor more design constraints, wherein at least one of the one or moredesign constraints comprise: a maximum slope constraint for a firsttriangle of the two or more triangles in the triangular surface mesh,wherein the maximum slope constraint comprises a maximum angle between anormal vector of the first triangle and a reference vector; (e) the CADapplication projecting, in the computer, one or more heights of thevertices of the first triangle onto the one or more design constraintsets such that the normal vector satisfies all of the one or more designconstraints, wherein the projecting comprises modifying the one or moreheights by a minimum Euclidian distance, and wherein the modifyingchanges the triangular surface mesh into a converged triangular surfacemesh; and (d) a design of the surface represented by the convergedtriangular surface mesh, wherein the design is generated and rendered inthe CAD application based on the projecting, wherein the surface isconstructed based on the design.
 10. The system of claim 9, wherein theone or more design constraints further comprise: a surface orientedminimum slope constraint for the first triangle, wherein the surfaceoriented minimum slope constraint comprises a minimum angle between thenormal vector of the first triangle and a directional vector.
 11. Thesystem of claim 9, wherein: the surface comprises a land surface; andthe land surface is constructed based on the design.
 12. The system ofclaim 11, wherein: the design comprises a computer-aided design (CAD)for the land surface.
 13. The system of claim 9, wherein the referencevector comprises a Z-vector representative of gravity.
 14. The system ofclaim 9, wherein: the application performs the projecting iterativelyuntil all of the one or more design constraints are satisfied.
 15. Thesystem of claim 9, wherein: the application projects onto the designconstraint set for the maximum slope constraint by moving the normalvector onto a cone defined by an origin of the reference vector on thetriangular surface mesh and a unit circle around the reference vectorthat is defined by the maximum angle.
 16. The system of claim 9,wherein: the application projects onto to the design constraint set forthe maximum slope constraint by moving the normal vector onto apolygonal shaped pyramid defined by the reference vector on thetriangular surface mesh and a regular polygon around the referencevector, wherein edges of the polygonal shaped pyramid are at the maximumangle with the reference vector.